© ESO 2018

1 Introduction

The possibility of determining for the first time the absolute proper motions of stars in the satellites of the Milky Way opens up a whole new window for understanding their dynamics, origin, and evolution, as well as that of the Milky Way itself. The data presented in the Second Gaia Data Release (hereafter DR2, Gaia Collaboration 2018b) allows us to achieve this goal. In this paper we study the proper motions (PM hereafter) of stars in a large sample of globular clusters, in the classical dwarf spheroidal galaxies and one ultra-faint system, and in the Large and Small Magellanic Clouds (LMC and SMC hereafter).

A plethora of interesting science questions can be addressed with this dataset. In this Introduction, we do not aim to be fully comprehensive, but we mention a few topics to set the context, to highlight the power of the unprecedentedly accurate absolute PM measurements, and also to fan curiosity in the community for exploring this outstanding dataset themselves.

Proper motion studies of satellite systems, such as the globular clusters and dwarf galaxies of the Milky Way, have a long history, starting from the use of photographic plates that were sometimes taken with a time baseline longer than 100 years (see Meylan & Heggie 1997 and van Leeuwen et al. 2000 for interesting and thorough historical reviews on the determination of PM of stars in globular clusters). More recently, the space missions HIPPARCOS and the Hubble Space Telescope (HST), and of course the Gaia mission in its first data release (Gaia Collaboration 2016), have demonstrated the enormous power of space-based astrometry. HIPPARCOS data (Perryman et al. 1997) have been used for many purposes, and in particular, for studying the dynamics of nearby open clusters (e.g. van Leeuwen 1999, 2009), and although HIPPARCOS did not observe stars in globular clusters, it provided an absolute reference frame that was used to derive the orbits of 15 globular clusters from photographic plates, for example (Odenkirchen et al. 1997). On the other hand, the HST has carried out several large (legacy) surveys (e.g. Soto et al. 2017) that have allowed studies of the dynamics of globular clusters and of the Milky Way satellites, and it has even constrained the motions of our largest neighbouring galaxy M31 (Sohn et al. 2012). In all these cases, relative astrometry is done using background quasars and distant galaxies to define a reference frame, and typically, a time baseline of 5–10 yr is used. This has been a highly successful approach, and has, forexample, allowed researchers to develop the idea that the Magellanic Clouds may be on their first infall (Kallivayalil et al. 2006b; Besla et al. 2007), to place constraints on the mass of the Milky Way from its most distant satellite Leo I (Boylan-Kolchin et al. 2013), and also to argue in support of the conjecture that dwarf galaxy satellites may lie on a vast polar plane based on the first constraints on their orbits (Pawlowski & Kroupa 2013).

This brief overview gives a flavour of the palette of scientific results that can be derived from accurate PM information of the satellites of the Milky Way. In combination with knowledge of the line-of-sight velocities, PM can be used to derive orbits for these systems. This is interesting for very many reasons, some of which we highlight below.

The orbits of globular clusters can shed light on their formation and evolution, for example, which may have formed in situ and which could be accreted (Searle & Zinn 1978; Mackey & Gilmore 2004; Renaud et al. 2017). Furthermore, knowledge of the orbits helps understanding the effect of tides and the interplay with internal processes, such as evaporation, mass segregation, and two-body relaxation (Djorgovski & Meylan 1994; Baumgardt & Makino 2003). Based on the orbits it is also possible to aid the search for extra-tidal stars and streamers, which are very useful for constraining the gravitational potential of the Milky Way because of the coldness of such streams (Küpper et al. 2015).

In the case of the dwarf galaxy satellites of the Milky Way, knowledge of the orbits also has multiple implications that range from the scale of the formation of the smallest galaxies in the Universe to constraints and challenges to the cosmological model. By determining the orbits of dwarf galaxies, we can establish the effect of the environment on their evolution, including star formation and chemical enrichment histories (Tolstoy et al. 2009), and also the effect of ram pressure stripping, and we can place constraints on the hot gaseous halo of the Milky Way (Nichols & Bland-Hawthorn 2011). The structure of these small galaxies may also have been strongly affected by tidal interactions with the Milky Way, and to quantify the importance of this process, knowledge of the orbits is imperative (Kazantzidis et al. 2011). Furthermore, such knowledge also allows to establish whether there is internal rotation and its amplitude (Battaglia et al. 2008), which is relevant for understanding the formation path of dwarf spheroidal (dSph) galaxies. For the ultra-faint galaxies, whose nature is debated, PM are also useful to identify interlopers, which is particularly important for establishing whether these systems are (on the verge of being) disrupted or embedded in a dark matter halo.

The orbits of the Milky Way satellites (both globular clusters and dwarf galaxies) also provide information on the Milky Way itself, such as its dynamical mass (e.g. Wilkinson & Evans 1999). It is likely that the internal dynamics of the Milky Way have also been affected by the gravitational influence of, in particular, the Sagittarius dwarf (Gómez et al. 2013) and the LMC (Bekki 2012; Gómez et al. 2015), and improved knowledge of the orbits of these objects will allow us to understand what their effect has been. On the other hand, orbits also allow us to gain insight into how a galaxy acquires its satellite population. For example, it has been argued that the satellites lie preferentially on streams (Lynden-Bell & Lynden-Bell 1995), on a thin plane (Kroupa et al. 2005), or that they have fallen in groups (Li & Helmi 2008), of which the LMC/SMC and their recently discovered satellites are direct proof (Bechtol et al. 2015; Koposov et al. 2015). The Gaia DR2 data will allow us to establish how real and important these associations are, and also whether the orbits found are consistent with the expectations from the concordance cosmological model.

In this paper we analyse 75 globular clusters in our Galaxy, and we demonstrate that the Gaia DR2 PM measurements for these clusters are of outstanding quality, with the formal and systematic uncertainties being effectively negligible. In comparison to previous efforts (e.g. Dinescu et al. 2003; Casetti-Dinescu et al. 2007, 2010, 2013), the errors are reduced by nearly two orders of magnitude. This dramatic improvement will also enable detailed studies of the internal dynamics that could shed light onto how these objects formed and their evolutionary path (Gratton et al. 2012). Some of the questions that might be addressed include whether globular clusters have formed in mini-halos or are fully devoid of dark matter (Ibata et al. 2013). Do they host intermediate mass black holes (Baumgardt 2017)? Are there dynamical differences between the different populations known to be present in many globular clusters (Bellazzini et al. 2012; Bellini et al. 2015; Vesperini et al. 2013)? Has the formation process and evolution for in situ clusters been the same as for those that have been accreted? Have these processes left an imprint on the internal phase-space distribution of their stars? How many clusters show rotation, and what is the link to how they have formed (Hénault-Brunet et al. 2015)? Many of the globular clusters are also being targeted by radial velocity surveys (e.g. Lardo et al. 2015; Kamann et al. 2018), and the combination of Gaia DR2 with such datasets will be extremely powerful.

We also study the Magellanic Clouds, the nine classical dSph, and include the ultra-faint dwarf (UFD) Bootes I as an example of what can be achieved with Gaia DR2 data. Even though the dwarf galaxies are on average farther away, their mean PMs can be very well determined using Gaia DR2, and they are still above the systematic level. Although for many objects, the uncertainties are comparable to those achievable using the HST, the advantage of having a full view of these galaxies and of the PMs being in an absolute reference frame cannot be over-emphasised. For the dSph, establishing their internal dynamics using this dataset is not yet feasible, however, although perhaps the combination of Gaia and the HST will allow to make progress before the end of the Gaia mission (as recently demonstrated by Massari et al. 2018). For the Magellanic Clouds, Gaia DR2 gives a clearer, more detailed view of the internal dynamics than has ever been possible before, with measured PMs for millions of sources.

The paper is structured as follows. The main part introduces the DR2 data, methods, and analysis, including orbit integrations, and details are given in the appendix. The appendix also contains tables with the measured PM for the objects we studied, as well as a list of the orbital parameters we derived. More specifically, in Sect. 2 of the main paper we present the Gaia DR2 data, with emphasis on the astrometry, the selection procedures, and the methods. Section 2.1 focuses on deriving the proper motions of the globular clusters and dSph, and in Sect. 2.2 we describe the procedures that are tailored for the LMC and SMC. We then present the various analyses of the datasets that we have carried out, and which allow us to show the superb quality of the data. Section 3 concentrates on the globular clusters, Sect. 4 on the dSph, and Sect. 5 on the Magellanic Clouds. In Sect. 6 we determine the orbits of the satellites using different Galactic potentials, a showcase of the fantastic possibilities that Gaia DR2 offers for studies of the dynamics and origin of the satellites of the Milky Way. In Sect. 7 we discuss our findings, provide an example of the use of DR2 astrometry to find tidal debris, present a summary of what lies beyond a straightforward analysis of the data such as that presented here, and also what will need to wait for later Gaia data releases (i.e. the limitations of the Gaia DR2 dataset). We present our conclusions in Sect. 8.

2 Data and methods

The data we used are the second Gaia data release as described in Gaia Collaboration (2018b). Further details on its validation may be found in Arenou et al. (2018). The procedures to derive the Gaia astrometric solution (also known as AGIS) are described in detail in Lindegren et al. (2016, 2018). We recall that the astrometric parameters are absolute in the sense that they do not rely on an external reference frame.

2.1 Globular clusters and dwarf galaxies

The sample of globular clusters analysed in this paper includes half of the whole population of globular clusters in the Milky Way. We focus mostly on the clusters that are located within a distance limit of 12 to 13 kpc to achieve a reasonable compromise on the number of stars with reliable astrometric solutions. It is important to bear in mind that the astrometric solutions for stars in areas of high stellar density, such as the cores of the clusters, are more likely to be disturbed by image blending and onboard image selection. This plays a significant role when observing more distant clusters and affects the fainter stars in particular (see e.g. Pancino et al. 2017). Our selection also takes into account the ability of distinguishing (in PM and parallax space) the cluster stars from those in the field, both as a function of distance from the cluster centre and of magnitude. Furthermore, clusters at low galactic latitude have also generally been avoided to escape confusion with field stars. The top panel of Fig. 1 shows a histogram of the distance distribution for the 75 globular clusters. The bottom panel exemplifies how the standard uncertainties on the cluster PM in declination vary as a function of the number of cluster members used¹.

As Fig. 2 shows, we also studied the classical dSph and one UFD galaxy, Bootes I. UFD galaxies are intrinsically very faint, as their name indicates, and this implies that there are very few stars on the red giant branch (RGB), and depending on the distance to the system, there may be even fewer because of the somewhat bright faint magnitude limit of Gaia (G = 21). Bootes I is the best UFD case for Gaia DR2, because its RGB is relatively well populated (at least in comparison with other UFDs), and it is relatively near (at 60 kpc, Belokurov et al. 2006). These conditions allow us to apply a homogeneous selection and analysis procedure to all the dwarfs in our sample, which we find highly desirable at this point. With external knowledge of radial velocity members, for instance, it might be possible to derive the PM for more UFDs, but the Bootes I case already illustrates the problems to be faced with Gaia DR2 data for this type of system.

The selection procedure, which is described in detail in Appendix A.1, starts with the extraction of data for each object from the GACS archive. The archive provides us with the astrometric parameters, their standard uncertainties and error correlations, the photometric data with standard uncertainties (flux values and fluxes converted into magnitudes), various statistics relating to the astrometric and photometric solutions, and radial velocities where available for our analysis. Depending on the nature of the object analysed, we set different magnitude limits. For the dwarf galaxies, we first considered stars with G < 21. For the globular clusters, the limit was generally set at G = 20, but in a fewcases, we took a brighter value to limit the contamination by field stars. This was necessary for clusters at low galactic latitude in particular.

The Gaia sky coverage can locally show strong variations that can affect the selection of members with good astrometric solutions (Gaia Collaboration 2018b; Arenou et al. 2018). In addition, for many of the globular clusters, thecentral core is often poorly resolved. These conditions are reflected in the standard uncertainties of the derived parameters, but are unlikely to cause a systematic bias in the results. The most strongly affected cluster ω Cen (see Fig. A.6) still shows good astrometric data for very many stars. In the case of the dSph galaxies, the most affected object is Sextans, as can be seen from Fig. 2. The inhomogeneous distribution of sources isrelated to the number of independent scans in the field of view towards the dwarf. To determine the astrometric parameters reliably, a sufficiently high number of truly independent scans is necessary. This is measured by the parameter visibility-periods-used, which has to reach a value greater than 5 for a five-parameter solution for an object (i.e. including the PMs and parallax) to be considered reliable (Lindegren et al. 2018), otherwise, only its position on the sky is determined. There are other instrumental effects that affect the astrometric parameters, and these are discussed elsewhere in the paper and in Lindegren et al. (2018).

The first step in our procedure to derive the motions of the satellites is to focus on an area of the sky, centred on the assumed centre of the object of interest and with an assumed maximum radius. For the dwarf galaxies, these radii were fixed at 2 deg, except for the Sagittarius dwarf, for which we took 3 deg (we also excluded stars within one tidal radius of its nuclear globular cluster M54). For the globular clusters, we interactively explored the data using the TOPCAT software (Taylor 2005), and then made a pre-selection of members based on the concentration of the PMs, followed by a cutoff in parallax, as well as on inspection of the colour-magnitude diagram (for more details, see Appendix A.1).

Because of their low stellar density contrast and the consequently higher number of contaminants (non-member stars) in the field of view, we applied additional selection criteria for the dSph galaxies in order to obtain a more robust estimate of the mean PMs. First, we only considered stars within 1.5× the tidal radius (r_t) of each dwarf (taken from Irwin & Hatzidimitriou 1995; Roderick et al. 2016, for Bootes) except for Sagittarius, where we considered all the stars in the 3 deg radius field of view. Then for all dSph, we also performed a cut in relative parallax error to remove foreground sources, as nearby stars will have relatively good parallaxes, especially in comparison to the stars in the dwarf galaxies. The relative error we used is 0 < σ_ϖ∕ϖ < 0.5 (which is equivalent to ϖ − 2σ_ϖ > 0), and corresponds to removing stars within roughly 5 kpc from the Sun. Finally, we used the distribution of sources in the colour-magnitude diagram (CMD) to isolate the giant branch (RGB and HB), as shown in Fig. 3 with the blue lines. In the case of the Sagittarius dwarf, we used a slightly different selection and focused on the reddest part of the RGB. The reason for this is the very large foreground, which overlaps substantially with the bluer portions of the Sagittarius RGB.

The astrometric solution to derive the PMs and parallaxes for the globular clusters and the dwarf galaxies follows the proceduresdescribed in van Leeuwen (2009); Gaia Collaboration (2017, see also Appendix A.1). A joint solution for the PM and parallax is obtained that takes into account the full error correlation matrix as evaluated for each contributing star: $$\mathcal{N}={\mathcal{N}}_a+{\mathcal{N}}_v+{\mathcal{N}}_d.$$(1)

The three main contributions to the noise matrix in Eq. (1) come from the astrometric solution ${\mathcal{N}}_a$, the estimated contributions from the internal velocity dispersion on the PM dispersion ${\mathcal{N}}_v$, and the dispersion of the parallaxes from the depth of the cluster ${\mathcal{N}}_d$, respectively.For the dwarf galaxies, the second and third of these contributions could be ignored, as even the brightest stars in these systems still have standard uncertainties on the astrometric parameters that are relatively large in comparison². Although for most globular clusters the velocity dispersion shows a clear gradient with respect to distance from the cluster centre (see Fig. A.5), we did not take it into account. This would have required a detailed investigation of the actual distribution of the PMs as a function of radial distance, which is beyond the scope of the present paper. The internal velocity dispersion as implemented is an average over the cluster.

The procedure we used to determine the astrometric parameters is iterative and requires a first guess for the parallax and PMs. For the globular clusters, this first guess was obtained using the TOPCAT software (Taylor 2005), as described above. While iterating, several diagnostics are produced, and in particular, we plot the surface density as a function of distance from the centre of the cluster. Such a diagram often shows that the maximum radius initially considered in the data extraction step can be extended farther out (i.e. the background density has not yet been reached). In that case, we retrieved more data from the GACS archive using an increased radius and the latest values found for the PM and parallax. We then repeated the procedure, now with the starting guesses being those given by the latest astrometric solution. This process was repeated until it was clear that the maximum radius had been reached.

The maximum radius for the cluster, that is, the distance from the centre within which we still detect cluster stars (3σ from the mean PM, where σ is the error on the PM derived using Eq. (1)) was compared to the tidal radii r_t extracted from Harris (1996) and its 2010 update (Harris 2010, hereafter Harris10). Figure 4 shows that for the majority of the clusters, this maximum radius is between 1∕2 and 2 times the published estimate of the tidal radius. Clusters for which the maximum radius was found to be much smaller than r_t are often affected by a high-density field star population, making the detection of cluster members problematic. We note that r_t has typically been estimated by fitting a King profile to the projected density distribution of stars, and thus does not necessarily nor always reflect the true extent of a cluster (see e.g. Küpper et al. 2010).

In the case of the globular clusters, the contamination by field stars was checked through the dispersion diagrams (see Fig. A.4 for two examples), in which the distribution of PM and parallax was plotted against the standard uncertainties of the measurements, and compared with the expected distributions that include all noise contributions. A contaminating source, such as the SMC for 47 Tuc (NGC 104), shows as an offset over-density in one or more of these charts, and in that case was removed by applying a 3σ filter to the residuals in all three observables, that is, relative PMs and parallax.

We also note that the parallax reference value used for the data extraction was the Gaia parallax for the cluster. This can differ from what is considered the best value for the cluster based on the distance from the literature (see Sect. 3 for more details).

In the case of the dwarf galaxies, the iterative procedures are similar, except that further iterations with the GACS archive are not necessary given our choices of initial field sizes. We thus worked only with the data extracted in the first step, as described earlier in this section. We have found, however, that we obtained more reliable mean PM using only stars brighter than a magnitude limit in the range 19.1 < G < 20. This is the faintest magnitude at which the mean value of the astrometric parameters becomes stable and where the effects of contaminating field stars and the large measurement uncertainties of very faint stars are minimised.

Fig. 1 Top: distance distribution of the 75 globular clusters included in the present study. Bottom: standard uncertainties on the PM in declination as a function of number of cluster members nMemb used in the solution. The diagonal line represents a fit to the relation $\sigma {\mu }_{\delta }=a/\sqrt{nMemb}$ where we find a = 0.3 [mas/yr]. A similar dependence on the number of members is found for the parallax uncertainty (with a = 0.15 [mas]) and for σμα* (where a= 0.25 [mas/yr]).

Fig. 2 Field-of-view towards the dSph galaxies (the nine classical and one UFD) in our sample. The stars shown correspond to membersaccording to the photometric selection (on the RGB and BHB) and the astrometric procedure (within 2σ from the mean PM of the object). The striping apparent in Sagittarius and Sextans is driven in part by the scanning law. The hole in the centre of Sagittarius corresponds to the location of the globular cluster NGC 6715 (M54).

Fig. 3 Colour-magnitude of the stars in the field of view towards the dSph galaxies (the eight classical and one UFD) in our sample. The blue lines mark our (relatively tight) pre-selection of tentative members (on the RGB and BHB) that is fed to the pipeline to derive mean PMs. The coloured points indicate stars within 3σ of our determination of the mean PM of the object. This means that cyan points satisfy both the PM and CMD selections.

Fig. 4 Comparison between the tidal radii rt (according to Harris10) of the 75 globular clusters in our sample and the maximum radii at which we have been able to detect cluster members in the present study. The diagonal line represents the one-to-one relation.

2.2 Magellanic Clouds

The LMC and SMC present a different analytical challenge to the analysis of dwarfs and globular clusters, because they are very extended onthe sky and contain two orders of magnitudes more Gaia sources than any of the dwarfs or clusters analysed.

To simplify our analysis and ensure that the quoted (and plotted) PMs are relatively easy to interpret in terms of internal velocities, itis particularly helpful to define an orthographic projection of the usual celestial coordinates and PMs: $$\begin{array}{ll}x\hfill & =\cos \delta \sin (\alpha -{\alpha }_C)\hfill \\ y\hfill & =\sin \delta \cos {\delta }_C-\cos \delta \sin {\delta }_C\cos (\alpha -{\alpha }_C)\hfill \\ {\mu }_x\hfill & ={\mu }_{\alpha *}\cos (\alpha -{\alpha }_C)-{\mu }_{\delta }\sin \delta \sin (\alpha -{\alpha }_C)\hfill \\ {\mu }_y\hfill & ={\mu }_{\alpha *}\sin {\delta }_C\sin (\alpha -{\alpha }_C)\hfill \\ \hfill & +{\mu }_{\delta }\left(\cos \delta \cos {\delta }_C+\sin \delta \sin {\delta }_C\cos (\alpha -{\alpha }_C)\right).\hfill \end{array}$$(2)

The centres of the coordinate systems are chosen to be the dynamical centre of the HI gas for the LMC and SMC, (α_C,LMC, δ_C,LMC) = (78. °77, −69. °01) and (α_C,SMC, δ_C,SMC) = (16. °26, −72. °42) (Luks & Rohlfs 1992; Kim et al. 1998; Stanimirović et al. 2004)³. Figure 5 shows the density of stars in this x, y-plane for the LMC and SMC, with these centres assumed.

If we approximate each cloud as a thin disc with some bulk motion that rotates about a point with celestial coordinates (α_C, δ_C) with a constant angular velocity ω and no other streaming motion, it can be shown (Appendix B) that these coordinates are, to first order, straightforwardly related to the parameters that describe the position and motion of the disc. These approximations are reasonable towards the centre of the LMC, and serve as a first approximation for the SMC.

It is convenient to define n to be the unit vector normal to the disc (such that rotation is positive about n), with z the unit vector from the observer to the reference centre (α_C, δ_C) at the reference epoch. We then have the mutually orthogonal unit vectors in the plane of the disc l = z × n∕|z ×n| and m = n × l. These have the property that l points in the direction of the receding node (the intersection of the disc with the tangent plane of the celestial sphere).

When we define v_x, v_y to be the centre-of-mass motion of the cloud in the x and y directions and v_z to be the same along the line of sight (divided by the distance to the cloud, to put it in the same units) then we have, to first order, $$\begin{array}{ll}\partial {\mu }_x/\partial x\hfill & \approx a{v}_x-v_z+a{l}_x{m}_z\omega \hfill \\ \partial {\mu }_x/\partial y\hfill & \approx b{v}_x-n_z\omega +b{l}_x{m}_z\omega \hfill \\ \partial {\mu }_y/\partial x\hfill & \approx a{v}_y+n_z\omega +a{l}_y{m}_z\omega \hfill \\ \partial {\mu }_y/\partial y\hfill & \approx b{v}_y-v_z+b{l}_y{m}_z\omega \hfill \end{array}$$(3)

where with inclination i (the angle between the line-of-sight direction to the cloud centre and the rotation axis of the disc, with i > 90° for retrograde motion)⁴, and Ω the position angle of the receding node, measured from y towards x, that is, from north towards east, we have the components of l, m, and n being $$\left[\begin{array}{ccc}{l}_x& m_x& n_x\\ l_y& m_y& n_y\\ l_z& m_z& n_z\end{array}\right]=\left[\begin{array}{ccc}\sin \Omega & -\cos i\cos \Omega & \sin i\cos \Omega \\ \cos \Omega & \cos i\sin \Omega & -\sin i\sin \Omega \\ 0& \sin i& \cos i\end{array}\right],$$(4)

and $$a=\tan i\cos \Omega ,b=-\tan i\sin \Omega .$$(5)

This means that simply by finding a linear fit to the PM as a function of position on the sky, yielding the bulk motion perpendicular to the line of sight, and four gradients, we have four equations for four (in principle) free parameters: v_z, i, Ω, and ω. The first, v_z, produces a perspective contraction (or expansion) as the clouds appear to shrink as they move away from us (or the opposite). The last three describe the orientation and rotation of the disc, which also leave a signature in the PMs.

In practice, neither cloud is flat or expected to have perfectly circular streaming motion. The assumption of a constant angular velocity is approximately valid in the central few degrees of the LMC, but this breaks down at larger radii. Nonetheless, these approximations allow us to draw tentative conclusions about the orientation and velocity curve of the Cloud from these gradients that are simple to measure.

We could take some of the four “free” parameters from other studies, but in practice, we only ever did this for v_z. For the LMC, we took the line-of-sight velocity from van der Marel et al. (2002, 262.2± 3.4 km s⁻¹), and the distance from Freedman et al. (2001, 50.1 ± 2.5 kpc), and for the SMC, we took the line-of-sight velocity from Harris & Zaritsky (2006, 145.6± 0.6 km s⁻¹), and the distance from Cioni et al. (2000b, 62.8 ± 2.4 kpc). This gives us v_z,LMC = 1.104 ± 0.057 mas yr⁻¹ and v_z,SMC = 0.489 ± 0.019 mas yr⁻¹.

To determine the PMs of the Clouds, we selected sources using the following procedure:

• 1.

  To create a filter, we initially selected stars with $\rho =\sqrt{x^2+y^2}<\sin r_{\text{sel}}$ (r_sel = 5° for the LMC, r_sel = 3° for the SMC) and ϖ∕σ_ϖ < 10 (to minimise foreground contamination). We also selected only stars with G < 19 in this step to ensure that the spread in PM due to uncertainties is small compared to the difference between the PM of the Cloud and of the bulk of the foreground.

• 2.

  We determined the median PM of this sample, and preliminarily filtered on PM by removing any source where μ_x or μ_y lies morethan four times the robust scatter estimate⁵ of that PM component from the median.

• 3.

  We determined the covariance matrix of μ_x, μ_y for these stars, σ, and usedthis to define a filter on PM, requiring that μ^Tσ⁻¹μ < 9.21 to correspond to a 99% confidence region.

• 4.

  We applied this filter in PM, along with that in ϖ, to all stars with G < 20 within 8 degrees of the assumed centre of LMC or SMC to define our complete sample.

We iterated this procedure twice, first using the expected μ_x, μ_y given the quoted μ_α*, μ_δ. This gave us a median parallax for the stars in the two Clouds: − 19μas for the LMC, and − 0.9 μas for the SMC (compared to the expected values of ~ 20 μas and ~ 16 μas, respectively). This is consistent with the offset and variation reported in other sections of this paper and in Arenou et al. (2018). We then repeated the procedure using the values of μ_α*, μ_δ implied by the data, conditional on the source parallax taking this median value (taking into account the quoted uncertainties and correlations). This procedure left us with 8 million sources in the LMC and 1.4 million in the SMC.

Fig. 5 Density distribution on the sky of the stars selected as members of the LMC (left) and SMC (right). Positions are shown in the (x, y) coordinatesdescribed in the text (Eq. (2)). In all figures that use this coordinate system, the x-axis has been inverted so that it corresponds to the usual inversion of right ascension.

3 Analysis: Globular clusters

As described earlier, we have analysed 75 globular clusters, for which the data are presented in Table C.1. For each cluster we have derived the PM and parallax, and where data were available, the radial velocity.

3.1 First analysis and comparisons

Figure 6 compares the parallaxes derived from the Gaia data to those from the cluster distances given in Harris10. There is a systematic difference of –0.029 mas (the Gaia parallaxes being smaller), originating largely from the Gaia data, and a calibration noise level around that relation of 0.025 mas (Arenou et al. 2018; Lindegren et al. 2018). A small contribution might also come from the values given by Harris10. However, we have made a provisional check on these distance estimates using the Gaia photometric data by superimposing the HR diagrams for all the clusters using the distances and reddening values as presented in Harris10 (see Gaia Collaboration 2018a). We found that all the clusters are neatly aligned for the critical elements (mainly the position of the blue horizontal branch). This indicates that, as a group, the distance moduli and colour corrections are confirmed to be in mutual agreement to better than 0.1 magnitude.

The standard uncertainties, which measure the precision rather than the accuracy, of the cluster-parallax determinations are smaller or very much smaller than the overall calibration noise level, as shown in the top panel of Fig. 7. The actual errors on these parallax determinations are therefore dominated by the overall Gaia calibration noise and offset in the parallax values. As discussed in depth in Lindegren et al. (2018), these systematic errors are also apparent in the parallax distribution of quasi-stellar objects (QSOs, which reveal the same offset), and as we show in Appendix A.2, also in the parallaxes of stars in the LMC (localised fluctuations) and other dSph, and are due to the basic angle variation and scanning law of Gaia. It is therefore expected that their amplitude will be significantly smaller in future Gaia data releases. For the time being, and because the parallax uncertainties derived photometrically are smaller, we use the distances as given by Harris10 in the analyses that follow.

The observed PMs are mostly about one to two orders of magnitude larger than the parallaxes, and thus the measurements are very robust and significant (see the bottom panel of Fig. 7). A comparison with a series of studies (Dinescu et al. 1999, 2003; Casetti-Dinescu et al. 2007, 2010, 2013) is shown in Fig. 8, and indicates overall good agreement, and most notably that the errors have been reduced by nearly two orders of magnitude. It remains somewhat uncertain, however, if the same calibration noise level can be assumed for the PMs as for the parallax (but see e.g. Sect. 4.1). Nonetheless, this systematic will be much smaller than the amplitude of the PMs themselves.

Radial velocities as measured by Gaia (Cropper et al. 2018) are available for 57 of the 75 clusters, although there were 3 or more cluster stars with measured radial velocities for only 46 clusters. While future Gaia data releases will contain radial velocities for more of these sources, this highlights a need for dedicated high-precision spectroscopy of these clusters to properly complement the Gaia astrometry. Figure 9 shows a comparison between ground-based (from Harris10) and Gaia radial velocity measurements, indicating a good relation for clusters for which enough stars have spectroscopic measurements (darker points). The relation between the number of stars and the standard uncertainty on the mean cluster velocity indicates an average internal velocity dispersion of the order of 4 km s⁻¹. This estimate of the intrinsic velocity dispersions is very similar to what is observed for the PMs.

Figure 10 shows the distribution on the sky of the globular clusters in our sample, where the arrows indicate the direction of motion and the colour-coding reflects the amplitude of the tangential velocities. These were derived using the PMs listed in Table C.1 and the distances from Harris10.

Fig. 6 Comparison between parallaxes as derived from the Gaia DR2 data and parallaxes derived from the cluster distances as given in Harris10.

Fig. 7 Top: parallax error against the Gaia parallax as determined from the Gaia data for 75 globular clusters. The black curves are the 1 and 3σ limits. The vertical dashed line corresponds to 3 times our estimate of the systematic error on the parallax. Bottom: PM errors against the PMs in right ascension (open circles) and declination (solid circles) for the clusters in our sample. The curves represent the value of the PM for 100σμ. The PM measurement has a significance lower than 10σ only for NGC 6453, for which ${\mu }_{\alpha *}/{\sigma }_{{\mu }_{\alpha *}}~4$.

3.2 Further results from the globular cluster astrometric data

The outstanding quality of the Gaia DR2 data together with the absolute reference frame (free of expansion and rotation) in which the PMs are presented has also allowed us to clearly detect rotation in 5 of the 75 globular clusters in our sample. For 3 of these clusters (NGC 104, NGC 5139, and NGC 7078), this was already known (Bianchini et al. 2013), but we have also detected rotation in NGC 5904 and NGC 6656 (see e.g. the left panel of Fig. 11). An indication of rotation can also be observed in NGC 5272, NGC 6752, and NGC 6809. Similarly, Gaia data allow measuring expansion and contraction in globular clusters. Forexample, NGC 3201 (Fig. 11, middle) shows very clear perspective contraction, which is due to its very high radial velocity and relatively large parallax. From this we may determine the parallax of this cluster in the same way as this used to be done for the nearby Hyades open cluster (see van Leeuwen 2009, and references therein). The Gaia data as presented here for the radial velocity and the PMs thus provide a cluster parallax of 0.221 ± 0.0086 mas, at about 2σ from the value of 0.204 mas given by Harris10. Finally, for NGC 6397, a cluster considered to have been subject to core collapse, we can still see a signal of the expanding halo (Fig. 11, right), clearly different from the expected very weak perspective contraction signal.

Furthermore, we find that our clusters have velocity dispersion profiles that decline with radius (Fig. A.5), and that several clusters show a slight increase in the outskirts, probably as the result of a halo of more loosely bound stars (as evidenced also by their spatial extent, see e.g. Olszewski et al. 2009; Carballo-Bello et al. 2012; Navin et al. 2016; Kuzma et al. 2018). This increase is found at a distance where contamination by field stars should not yet be important.

Fig. 8 Comparison of the Gaia PMs (in right ascension: top, and declination: bottom) to measurements reported in Dinescu et al. (1999, 2003); Casetti-Dinescu et al. (2007, 2010, 2013) for 31 globular clusters.

4 Analysis: Dwarf spheroidal galaxies

The procedures described in Sect. 2.1 allow us to determine the mean PMs of the dSph in our sample. As discussed earlier, we focus in this paper on the classical dSph, and have included in our sample one example of an ultra-faint galaxy, Bootes I. The resulting mean ${\mu }_{\alpha }^{*}$ and μ_δ, as well as the astrometric parameters and uncertainties, are listed in Table C.2.

The efficiency of our selection procedure in removing most of the foreground contamination becomes clear in Fig. 12. Stars surviving the proposed criteria are shown as cyan dots for each of the dwarfs, and they clearly clump much more strongly in the diagrams than the likely non-members (shown as black points). In this figure the blue ellipses indicate the contours corresponding to 3σ dispersion around the mean ${\mu }_{\alpha }^{*}$ and μ_δ (σ is computed taking into account the covariances, using the standard error on the mean $\times \sqrt{N_{*}}$, where N_* is the number of stars used to measure the mean PMs).

The possibility of selecting members via their PM that the Gaia DR2 data provide opens a new window for understanding the structure and extent of the dSph. In particular, Fig. 2 shows that some of the dwarfs in our sample present spatial asymmetries (e.g. Fornax in the top right corner, the Sculptor outskirts appear somewhat boxy), while there is an indication of tidal streams in the case of the Carina dSph. A more detailed analysis of these features is beyond the scope of this paper, but the quality of the dataset certainly makes this possible.

Fig. 9 Top: comparison between the ground-based and Gaia radial velocities for 52 clusters with at least two measurements. Bottom: standard uncertainties on the mean radial velocities as a function of the number of stars contributing to the mean. The diagonal line shows the effect of an additional contribution of 4 km s−1 originating from the internal velocity dispersion.

4.1 Systematics, correlations, and dispersion

Table C.2 shows that the average parallax is negative for several dSph⁶. This systematic error, similar to that found for the globular clusters and the Magellanic Clouds (and the QSOs Lindegren et al. 2018), is present in the different fields, and its amplitude varies from object to object (for more details, see Appendix A.2). The average offset, computed as the difference between the expected parallax (based on the distances from McConnachie 2012), and the parallax from DR2 for all dSph, is − 0.056 mas⁷, and when Leo I is excluded, it is − 0.038 mas (the parallax offset is − 0.21 mas for Leo I).

The PM maps shown in Fig. 12 reveal that the dwarfs are extended in PM space. The main contributor of this dispersion is not intrinsic but is due to the uncertainties, which are typically very large. For example, for an object such as Sculptor, the individual velocity errors for G ~ 18 mag stars are of the order of 80 km s⁻¹ (for G ~ 20, they are > 200 km s⁻¹), compared to the expected internal dispersion of order of 10 km s⁻¹. Therefore, measuring the intrinsic dispersion for these systems does not appear to be feasible with the data provided by DR2. However, it may begin to become feasible with later data releases, and certainly with an extension of the Gaia mission.

The measurements of the PM are also affected by the scans and varying astrometric incompleteness, which introduce a pattern in the parallax and PM field that is only readily apparent for sufficiently large objects on the sky (see also Arenou et al. 2018; Lindegren et al. 2018). This is illustrated for the LMC and SMC in Figs. A.9 and 16, and it is also present for example for the Sagittarius dSph, as shown in Figs. 2 and 13.

For sufficiently large objects on the sky, the banding pattern is averaged out, and the mean PM is more robust⁸. The global dispersion is larger than expected just from random errors, however (because of the offset from bin to bin in the pattern). We quantify this effect in Fig. 14, where we show the distribution of the mean PM in α and δ computed in bins of 0.2 × 0.2 deg² size that contain at least 100 stars⁹, and after 2.5σ clipping to remove outliers. The fact that the mean value changes from bin to bin is at least partly caused by the finite number of stars in each bin, as well as by the random errors. We tried to estimate the residual systematic error by modelling this distribution assuming the mean PM derived using all member stars, and assuming that the errors are Gaussian. We drew a new PM from this mean for each star, assuming its quoted uncertainty and correlations, and recomputed the mean using all the stars in the bin (we also assumed an intrinsic velocity dispersion of 10 km s⁻¹). This is the red histogram in Fig. 14. Clearly, the observed distribution is wider (blue), and the difference between the two can be used to compute the systematic error as ${\sigma }_{mock}^2+{\sigma }_{sys}^2={\sigma }_{tot}^2$. We find this to be σ_sys ~ 0.030 mas yr⁻¹ and ~ 0.036 mas yr⁻¹ in right ascension and declination directions, respectively. These values are of slightly smaller amplitude than those derived by Lindegren et al. (2018) from a sample of QSOs.

For systems that are smaller on the sky, and in particular for those that would fall in a single bin, their PM may be offset by this much. Objects such as the larger dSph Fornax, Sculptor, etc., are likely not affected by this systematic (because it averages out), but for systems such as Leo II, it should be considered. In our subsequent analyses we thus considered the amplitude of the systematic uncertainty to be 0.035 mas yr⁻¹. Table C.2 shows that in many cases, this systematic error is larger than the random error on the measurement of the mean PM of a dSph.

Table C.2 also shows strong correlations in the different mean astrometric parameters derived for the dSph in our sample. These correlations vary from object to object in amplitude and direction (see Appendix A.2), and it is important to take them into account in the derivation of the orbital parameters, for instance.

Fig. 10 Sky distribution of the 75 globular clusters in our sample in Galactic coordinates. Their tangential velocities are denoted by the size and direction of the arrows, while the colours indicate their line-of-sight velocities. The inset shows a zoom-in of the central 60 × 60 deg2.

Fig. 11 PM “systematics” in radial (blue triangles) and transverse (red squares) directions, as a function of distance to the cluster centre. The cyan line is a fit to the variation of the mean radial component of the PMs as a function of distance, while the grey line is the expected trend resulting from perspective contraction or expansion. From left to right: NGC 6656 shows a strong rotation signal as well as perspective expansion; NGC 3201 shows no rotation but very strong perspective contraction; and NGC 6397 shows halo expansion associated with core collapse, clearly different from the signal expected from the small perspective contraction. The vertical line indicates the tidal radius of the cluster reported in Harris10.

4.2 Comparison to the literature

We have compared the PM we derived with our selection criteria and those we would obtain if we were to use only stars identified as members from publicly available radial velocity catalogues (from Armandroff et al. 1995; Kleyna et al. 2002; Muñoz et al. 2006; Battaglia et al. 2006, 2011; Sohn et al. 2007; Mateo et al. 2008; Walker et al. 2009, 2015). We have found very good agreement (i.e. the estimates differ by less than 1σ, see Arenou et al. 2018). The main disadvantage of using external information is its heterogeneous nature. Furthermore, this information is not available for all the dSph in our sample, and the sample of stars with radial velocities for any given dSph is typically smaller by a factor ~2 than is found with our selection and analysis, which are based exclusively on Gaia DR2 data, even down to the same magnitude limit.

Figure 15 compares our measurements of the PM of the dSph to astrometrically derived values reported in the literature. These are from Piatek et al. (2003, 2004) for Carina; from Pryor et al. (2015); Sohn et al. (2017) for Draco; from Piatek et al. (2005) for Ursa Minor; from Piatek et al. (2006); Sohn et al. (2017); Massari et al. (2018) for Sculptor; from Piatek et al. (2007) for Fornax; from Lépine et al. (2011); Piatek et al. (2016) for Leo II; from Sohn et al. (2013) for Leo I, and from Casetti-Dinescu et al. (2018) for Sextans. Most of these measurements were obtained from space, except for Sextans, whose measurements are based on Subaru imaging. We excluded measurements for Sagittarius and Bootes I for visualisation purposes (because their PMs are much larger than for the other, more distant dSph).

A striking difference between previous estimates of the dSph PM and those obtained using Gaia DR2 data is the extent of the error bars, particularly for the objects for which more than 400 (and up to several thousand) astrometric members have been identified, such as Carina, Ursa Minor, Fornax, Sculptor, and Draco. In many cases, our measured PMs are consistent with the literature values at the 2σ level (given the large error bars of the latter). For the most recent astrometric measurements with the HST (which are therefore typically based on a larger baseline), the values appear to be closer and consistent with each other (e.g. Leo I), especially when the systematic uncertainties are taken into account (e.g. Sculptor and Draco). In the case of Sagittarius, we find that the Gaia DR2 PM is consistent with that of Massari et al. (2013), although it is now much more accurate. We here present the first measurement of the PM of the UFD Bootes I.

Fig. 12 PMs of the stars in the field of view towards the different dSph galaxies in our sample. Members, defined as stars within 3σ of the mean measured PM and located in the expected region of the CMD, are shown in cyan. The green points correspond to stars that also fall in the CMD-selected box, but are not within 3σ of the systemic PM. Especially for Sculptor and Fornax, it is quite clear that there may be more members, but very likely, the large errors on the PMs of individual stars place them beyond the 3σ ellipse.

Fig. 13 Counts in 0.2 deg wide bins on the sky with at least 100 stars for stars in the Sagittarius dwarf (left) and the average PM in μα* (middle) and μδ (right) for each of these bins.

Fig. 14 Blue histograms: distribution of average PM in the bins shown in Fig. 13; red histograms: a model with the derived average PM convolved with the random errors provided by the astrometric solution. The difference shows the amplitude of the systematic on the mean value of the PM for a bin of this size.

5 Analysis: LMC and SMC

5.1 Basic analysis

In Fig. 16 we show the median PM, in the x, y coordinate system defined in Sect. 2.2, of sources that meet our membership criteria for the LMC and SMC, binned by position on the sky. This is a demonstration of the extraordinary precision of the Gaia PM measurements (see also Fig. 24). The rotation signature is clearly visible in the LMC, and trends in the PMs of stars in the SMC are visible as well. Figure 16 also serves as a demonstration of the shortcomings of this data release. The banding or striping of the PMs that we discussed above, which is associated with different scans and has been investigated by Lindegren et al. (2018), is also clearly visible. The parallaxes show this as well (Fig. A.9).

We can characterise the trends seen in Fig. 16, to first order, by the central values and gradients. We calculated these either directly as a least-squares fit to the data or as a least-squares fit to the median PMs calculated in 0. °04-by-0. °04 bins in x,y plane. This latter fit was performed to reflect the fact that the most important errors in this analysis are systematic and depend on position on the sky, but it becomes less appropriate towards larger radii as Poisson noise becomes more important. The differences between the values derived using these two methods give a sense of the scale of the uncertainty associated with the position-dependent systematic errors.

In Tables B.1 and B.2, we show the central values and gradients (Eq. (3)) for the PMs. We provide values for all sources within various angular radii (i.e. $\rho ={\sin }^{-1}{(x^2+y^2)}^{1/2}<{\rho }_{\text{max}}$) and for annuli. For both the LMC and SMC, we show the values of i, Ω, ω that we found when we performed a least-squares fit to the derived gradients under the assumption that v_z takes the value implied by the known values of the line-of-sight velocity and distance to the Clouds. The results are reasonably consistent with one another (allowing for the fact that we expect ω to decrease farther out).

For the LMC we also give the implied values of v_z, i, Ω, ω when we place no constraint on the line-of-sight velocity. The value of v_z that this implies is of the order of 1.3– 2.0 mas yr⁻¹, which is similar to (but somewhat higher than) the value 1.104 ± 0.057 mas yr⁻¹ that is expected given the measured line-of-sight velocity and distance of the LMC. The effect of the line-of-sight velocity of the Clouds is to produce a perspective shrinking of the Cloud on the sky (a negative contribution to both ∂μ_x ∕∂x and ∂μ_y ∕∂y in Eq. (3)), similar to the effect seen in NGC6656 or NGC3201 (Fig. 11). The mismatch between the value derived from Gaia astrometry and that found from spectroscopy might be related to an actual contraction of the LMC disc (similar to the apparent expansion of NGC6397 in Fig. 11). However, the orientation of the LMC disc also plays an important role in these values, and the values of i derived when v_z for the LMC is fixed lie closer to those found in photometric studies, which tend to be in the range 25− 40°.

For the SMC, directly inverting Eq. (3) gives line-of-sight velocities that are completely inconsistent with those measured from spectroscopy (~ −0.8 mas yr⁻¹ as opposed to 0.489 ± 0.019 mas yr⁻¹). This may be due to the inadequacy of modelling the SMC as a flat disc, or a real expansion of the SMC (which, again, is degenerate with line-of-sight motion). However, forcing v_z to take the value expected from the measured distance to the SMC and its line-of-sight velocity gives us a model that has a disc inclination ~74°, which is broadly similar to that measured for the Cepheid population (64.4° ± 0.7°: Subramanian & Subramaniam 2015).

In Fig. 17 we show the residual PMs after we subtracted off a gradient in PM corresponding to our first-order approximation, with the parameters v_z, i, Ω, ω found for sources within angular radii ρ_max = 3° of the centre for the LMC, and ρ_max = 2° of the centre for the SMC. This shows the scale of the spatially correlated errors in PM more clearly, and in Fig. A.10 we show the variation in 1D stripes across the LMC (as well as the variation of the parallaxes), to allow an easier quantification. The residuals are comparable to those found in Sect. 4.1. The residuals in the centre are rather small, but become larger far from the centre for the LMC, in the opposite sense to the variation from the median shown in Fig. 16. This is because of our assumption of constant ω, which breaks down badly at large radii, as the rotation curve becomes flat.

Figure 17 shows indications of the impact of the LMC bar on the kinematics of the disc. The residual PM near the upper (as we see it) side of the bar tend to be negative, while those on the lower side tend to be positive. This indicates that the stars aremoving at faster-than-circular velocities at these points; this is consistent with stars belonging to the x1 orbit family, which is elongated along the bar (Contopoulos & Papayannopoulos 1980).

Fig. 15 Comparison to literature values of the proper motions of the dSph. Symbols with the same colour correspond to the same dSph, where filled symbols with error bars are those derived in this paper, and open symbols surrounded by ellipses correspond to the literature values. The error bars have the size of ${ϵ}_{{\mu }_{\alpha }*}$ and ${ϵ}_{{\mu }_{\delta }}$, as reportedin Table C.2. The black cross in the bottom left corner indicates our estimate of the systematic uncertainty on the PMs.

Fig. 16 PMs of stars in the LMC (upper) and SMC (lower) showing the components μx (left) and μy (right), which are described in the text (Eq. (3)). The colour shows the median PM in each pixel (after filtering). The black density contours are logarithmically spaced, such that the outermost contour is at a source density 100 times lower than the highest density. The centre of each colour bar is chosen to be the median PM of all sources. This is not the same as the PM derived for the clouds below, because the sources are not distributed symmetrically around the assumed dynamical centre; the photometric and dynamical (from the HI disc) centres are offset from one another. Trends in PM, particularly the trend associated with rotation in the LMC, are clearly visible. The banding associated with the Gaia scanning law, and as seen in the parallaxes, are clearly visible as well.

5.2 Uncertainties and comparison to the literature

The uncertainties on the measurement of the centre-of-mass motion of the LMC and SMC using Gaia data are completely dominated by systematic, rather than random, uncertainties. The estimates of these quantities that we show in Tables B.1 and B.2 are consistent to around the 10 μas level. This is smaller than the systematic uncertainty on PMs calculated in Sect. 4 or that on a large scale (of a few tens of degrees) derived from the PMs of quasars observed by Gaia (Lindegren et al. 2018), which is ~ 28 μas yr⁻¹ in each component. There is no clear choice of the correct values to take from Tables B.1 and B.2. However, we can note that the approximation of constant angular velocity becomes poor for the LMC beyond about 3°, and that the density of stars in the SMC is small beyond a similar radius, so this appears to be a sensible choice. We therefore adopted the values we found using all stars within these radii as our best estimates. These are (μ_α*,0,LMC, μ_δ,0,LMC) = (1.850 ± 0.030, 0.234 ± 0.030) mas yr⁻¹ and (μ_α*,0,SMC, μ_δ,0,SMC) = (0.797 ± 0.030, −1.220 ± 0.030) mas yr⁻¹, where our uncertainty was estimated from the ~10 μas yr⁻¹ variation listed in Tables B.1 and B.2 and the ~ 28 μas yr⁻¹ large-scale systematic uncertainty.

Kallivayalil et al. (2013) give an overview of recent estimates of the PMs of the LMC and SMC, including their own, found using HST three-epoch astrometry. Our estimates are consistent with theirs and with almost all of the values they cite, as well as with the values found by van der Marel & Sahlmann (2016) using PMs found from the Tycho-Gaia Astrometric Solution (Gaia Collaboration 2016; Lindegren et al. 2016), which were (1.872 ± 0.045, 0.224 ± 0.054) mas yr⁻¹ and (0.874 ± 0.066, − 1.229 ± 0.047) mas yr⁻¹ for the LMC and SMC, respectively.

It is worth noting that the measured centre-of-mass PM is dependent on the chosen (or derived) centre of the Clouds. The centre of the HI gas disc of the LMC, which we assume to be the dynamical centre of the Cloud, is close to, but not exactly the same as, the centre that was derived from HST PMs by van der Marel & Kallivayalil (2014), which was that assumed by Kallivayalil et al. (2013) and van der Marel & Sahlmann (2016).

Instead of assuming that the stars’ dynamical centres lie at the dynamical centre of the HI gas disc in each case, we could assume that they lie at the photometric centres of the Clouds at (α_C,LMC,phot, δ_C,LMC,phot) = (81. °28, −69. °78) (van der Marel 2001) and (α_C,SMC,phot, δ_C,SMC,phot) = (12. °80, −73. °15) (Cioni et al. 2000a), respectively. When we do so, we derive mean PMs (μ_{α*,0,phot, LMC}, μ_{δ,0,phot, LMC}) = (1.890, 0.314) mas yr⁻¹ and (μ_{α*,0,phot, SMC}, μ_{δ,0,phot, SMC}) = (0.685, −1.230) mas yr⁻¹. Therefore, this variation is stronger than the variation due to the large-scale systematic uncertainties.

Our uncertainties, including systematics, are comparable with those of previous studies. The sheer number of sources spread across the Clouds for which Gaia provides accurate PMs is extraordinary, and this will allow astronomers using DR2 to make a detailed mapping of the dynamics of the Magellanic Clouds.

Fig. 17 Residual PM, after subtraction of a model PM field for the LMC (top) and the SMC (bottom). For both objects, the majority of the variation disappears, and the banding in PM is more clearly visible. Black density contours are spaced in the same way as in Fig. 16. For the LMC, the imprint of orbits on the bar can be seen as the bluer area on the lower side of the bar (as it appears in the plot) and the redder areaon the upper side of the bar. The model that is subtracted is fit from sources within angular radii ρmax = 3° of the centre for the LMC, and ρmax = 2° of the centre for the SMC.

Fig. 18 Rotation curve (top) and median vR (bottom) of the LMC. The assumed values for the centre-of-mass velocity and orientation of the disc (i and Ω) are taken from a fit to all stars within angular radii ρ < 3° of the LMC centre. Angular distances and velocities given on the lower and left axes have been converted to real-space values on the upper and right axes assuming a distance to the LMC of 50.1 kpc (Freedman et al. 2001). The points shown in the upper panel are derived from observed line-of-sight velocities of old and young stars by van der Marel & Kallivayalil (2014, their Table 4.)

5.3 Rotation curve

If we assume that we know the orientation of the LMC (or SMC), and that the motions we see are confined to the plane, it is possible to de-project the observed motions (minus the bulk motion) onto that plane. In Appendix B we give the mathematical details of how this was performed. The SMC is less suitable for approximation as a simple flat rotating disc than the LMC, therefore we did not attempt this here. Tentative evidence that there is some sense of rotation of the SMC stars is provided by the consistent, but small, measurement of ω for different annuli of stars shown in Table B.2.

In Fig. 18 we show the resulting median tangential velocity, v_T (the rotation curve) and median v_R as a function of de-projected radius R for the LMC(note that R, v_R, and v_T are de-projected position and velocity, i.e. in the plane of the LMC), with v_x, v_y, i, Ω as determined from a least-squares fit to the filtered data for angular radii ρ < 3°, holding v_z fixed. The figures show the median value of v_T (or v_R) as we increase R, in non-overlapping bins of 40 000 sources. In both cases we also divide the sample into sources with y > 0 and y < 0 as a consistency check. We also show the velocity curve derived from line-of-sight velocities by van der Marel & Kallivayalil (2014), separately for young and old stars. This shows that the precision of the PM rotation curve from Gaia is competitive with those in line-of-sight velocity curves (which are derived from very many more sources). The Gaia data contain both old and young stars, therefore it is expected that the Gaia rotation curve lies between the two curves from the old and the young populations.

The rotation curve for the LMC rises approximately linearly for R ≲ 3° (which provides post hoc motivation for us to choose stars within a projected angular radius ρ of 3° to determine the other parameters of the disc; the assumption of constant ω is reasonable over this radius). The dense coverage of the LMC provided by Gaia allows us to resolve this rise in the rotation curve in a way that was not possible with the relatively sparse coverage provided by the HST (van der Marel & Kallivayalil 2014). The scale over which the rise occurs is closer to that found for the old stellar population using line-of-sight velocities by van der Marel & Kallivayalil (2014). We also show the rotation curve that we derived when we allowed v_z to vary freely, which is broadly similar.

As a further sanity check, we divided the sources into blue and red bins, with the division at G_BP − G_RP = 0.6, and we plot the rotation curve in each case. As one would expect, the blue stars tend to have a higher rotation velocity out to about 4°, by about 5−10 km s⁻¹, reflecting a lower asymmetric drift.

The v_R plot is rather harder to interpret. It is pleasingly close to 0 km s⁻¹ in the inner regions, and the variation we see at projected radii < 2° may well be related to the effect of the bar. In the outer regions there is a significant difference between the trends seen at positive and negative y-values, which is notseen in the rotation curve. This might be due to non-equilibrium effects that are possibly caused by the past interaction of the LMC and SMC (e.g. Besla et al. 2016), but further interpretation is beyond the scope of this paper.

6 Orbital Integrations

In this section we present the results of the orbital integrations for the globular clusters in our sample and for the dSph. We did not perform integrations for the Magellanic Clouds as this would require consideration of dynamical friction, which introduces additional degrees of freedom such as the total mass of each Cloud (e.g. Kallivayalil et al. 2013). For each object we integrated an ensemble of 1000 orbits whose initial conditions were drawn via Monte Carlo sampling the measurements and their uncertainties given in Tables C.1 and C.2 (i.e. we used the full covariance matrix and assume Gaussian errors). We produced two sets of Monte Carlo samples, considering only the random error on the observables, and considering in addition a systematic error of amplitude 0.035 mas yr⁻¹ (as estimated in Sect. 4.1) for each PM component. We then transformed these coordinates into Cartesian positions and velocities. Because the Gaia DR2 parallaxes (or distances) have larger (systematic) uncertainties than the measurements available in the literature (see Sects. 3 and 4), we used published values from Harris10 and McConnachie (2012) for the globular clusters and the dSph, respectively. The radial velocities were also taken from these databases. For the distance errors, we assumed an uncertainty in the distance modulus of ~ 0.05 mag, which corresponds to a relative distance error of ~0.023.

The initial conditions for the orbit integrations are listed in Tables C.3 and C.4 for the globular clusters and the dwarfs, respectively. In these tables we give the uncertainty of the positions and velocities derived from the 16th and and 84th quantiles (which would be a 1σ deviation if the distributions were Gaussian), which were obtained by marginalising over the other coordinates. However, as started earlier, we took into account the full covariance matrix in the orbit integrations. We considered three different Galactic potentials (labelled Model-1, -2, and -3) for the orbital integrations, and each of these potentials corresponds to a model previously published in the literature. Our goal was to understand how different their predictions are, and also how robust the conclusions. In the near future, it will be possible to use the data to constrain the model parameters and to understand which model performs best, for example by imposing self-consistency, but this exercise is beyond the scope of this paper.

Because each potential has its own set of characteristic parameters whose values were derived by fitting different observables, they each assume different values for the position of the Sun and its peculiar velocity as well as for the motion of the local standard of rest. This implies for example that the angular momentum of the dwarf galaxies and globular clusters can be slightly different for the various potentials.

6.1 Description of the Galactic potentials

The gravitational potential in Model-1 is axisymmetric, consisting of a stellar bulge and discs (which have a combined stellar mass of 5.4 × 10¹⁰ M_⊙), two gas discs (with a combined mass of 1.2 × 10¹⁰ M_⊙), and a Navarro et al. (1996, NFW) dark matter halo. The virial mass is 1.37× 10¹² M_⊙. The model was found using a Bayesian analysis of kinematic tracers by McMillan (2017, cf. their Table 3).

The Model-2 gravitational potential is axisymmetric, consists of a stellar bulge and disc modelled as Miyamoto-Nagai potentials (which have a combined stellar mass of 7.55× 10¹⁰ M_⊙), and a spherical dark matter halo.The mass at R < 200 kpc is 1.9 × 10¹² M_⊙. The model is the one of Allen & Santillan (1991) but with revised parameters from Irrgang et al. (2013, their Table 1).

The third potential we considered (Model-3) is the non-axisymmetric mass model described in Robin et al. (2003, 2012). The properties of the axisymmetric components (thin and thick stellar discs, stellar halo, interstellar matter disc, dark matter halo) are described in Robin et al. (2003), while those of the rotating non-axisymmetric bar, renormalised to have a total mass of 6.7 × 10⁹ M_⊙, are described in Robin et al. (2012). The mass of the Galaxy for R < 100 kpc is ~1.2 × 10¹² M_⊙. The computation of the potential is described in Bienayme et al. (1987) for the axisymmetric components and in Fernandez-Trincado et al. (in prep.) for the bar and halo components.

The three potentials have very similar mass distributions between ~ 3 and 40 kpc as measured by the circular velocity curves (shown in Fig. D.1). However, the potentials differ substantially both in the inner and in the outer regions. This will lead to some differences in the orbits and their characteristics, as we show below.

6.2 Results for the globular clusters

Fig. 19 shows some examples of orbits for the globular clusters integrated backward in time for 0.25 Gyr. For all the clusters, the orbits are very similar initially (at least for one orbital period), but they then begin to diverge, reflecting the differences in the Galactic potentials that were used (see also Fig. D.1). The differences for clusters that penetrate the regions dominated by the bar are particularly large; this is modelled as a non-axisymmetric component in Model-3 (in red), but not in Models 1 and 2 (in blue and black, respectively). Furthermore, some clusters appear to be on resonant orbits, and interestingly, as shown in the second panel of Fig. 19 for NGC 6441, this is true for the three potentials. A quick exploration reveals that several other clusters appear to be on similar types of resonant orbits. Some clusters, on the other hand, seem to be on chaotic orbits, but further analysis is required to establish this reliably.

The distribution of some of the orbital parameters for the globular clusters in our sample is shown in Fig. 20 and summarised in Table D.1. In these plots, the solid circles show the median value of the time averages over 10 Gyr of integration, while the error bars correspond to the 16th and 84th percentiles derived using only the statistical (and not the systematic) errors on the observables. These figures show that the orbits of globular clusters in our sample are very centrally concentrated, as most have their orbital apocentres within 10 kpc. Most of our clusters are on prograde orbits, and a small fraction are retrograde. We recall that this sample is focused on the inner 20 kpc, so we cannot establish with this dataset whether this is also the case for the outer halo clusters. The clusters with the highest L_z are NGC 3201 (retrograde) and NGC 4590 (prograde) are not shown in this figure, nor is NCG 5466 because of its large apocentre (~ 50 kpc).

In Fig. 20 the eccentricity is defined as the time-average of (r_apo − r_peri)∕(r_apo + r_peri). The eccentricity distribution is extended and rather uniform with many clusters on relatively radial orbits. The apparent trend that clusters with larger apocentres are on more radial orbits may be due to our sample selection (currently located within ~ 20 kpc from the Sun). A cluster with a large apocentre will almost only be included in our sample if it has a relatively radial orbit. It will be possible to draw more reliable conclusions about the distribution of orbital parameters when the whole globular cluster population has been analysed. Nonetheless, Fig. 20 already highlights that the greater differences in the orbital properties for the various potentials arise for clusters with apocentres in approximately the inner 5 kpc. These differences are in manycases larger than the error bars. This implies that there is room for improvement by performing a self-consistent dynamical model of the globular cluster population and the mass distribution in our Galaxy (e.g. Binney & Wong 2017).

Fig. 19 Examples of the orbits of some of the globular clusters in our sample. The different colours correspond to the different potentials: Model-1 (based on McMillan 2017, in blue), Model-2 (based on Allen & Santillan 1991, in black), and Model-3 (based on Robin et al. 2003, 2012, in red). The orbits of clusters that remain in the inner few kpc are quite different for the various potentials, while as expected, the differences are much smaller for those that have pericentres greater than ~ 2 kpc (e.g. NGC6496, fourth row). In these cases, the location of streams, if present, can be predicted much more reliably.

Fig. 20 Distribution of orbital parameters for the globular clusters in our sample. Three globular clusters are not shown in these plots: NGC 3201 and NGC 4590, which have apocentres ~ 30 kpc and large retrograde and prograde motions, respectively (|Lz| > 2000 kpc km s−1), and NGC 5466, which has an apocentre ~50 kpc (given its orbital parameters, it was likely associated with the Sagittarius dwarf, as proposed by Bellazzini et al. 2003). We note the concentration of clusters with small pericentre and apocentre around Lz ~ 0 for Model-3. Such clusters are on box orbits in the barred potential of Model-3, and hence do not conserve Lz, whose time average over the 10 Gyr of integration is approximately zero. The error bars correspond to the 16th and 84th quantiles obtained from the orbit integration using the statistical errors. The effect of the systematic error on the orbits of the globular clusters is found to be negligible and is not shown here.

6.3 Results for the dwarf galaxies

In Fig. 21 we plot the orbits of the different dwarf (spheroidal) galaxies. This figure reveals for example that the orbits of Draco and Ursa Minor look similar, and that the orbital planes of most dSph are different, with the orbit of Sagittarius, for instance, being orthogonal to those of Draco and Ursa Minor.

Figure 22 shows the distribution of orbital parameters for the dSph (and listed in Table D.2). We have plotted here the median values (over 10 Gyr of integration) and the 16th and 84th percentile range as the symbols with solid error bars (derived from the 1000 Monte Carlo realisations of the observables and their uncertainties). Using the same symbols, but now with dotted lines, we have plotted the orbital parameters derived including the effect of the 0.035 mas yr⁻¹ systematic uncertainty on each of the PM components. In general, the effects of the systematic errors on the characteristic properties of the orbits are relatively small. We note, however, that in both cases these error bars do not always properly reflect the uncertainties on the orbital parameters because of degeneracies. We include examples of the Monte Carlo realisations in the appendix (Fig. D.3) to give examples of these degeneracies.

Figure 22 reveals that most satellites are on (slightly) prograde orbits, while Fornax is retrograde (and possibly Leo II as well, although it is also consistent with prograde at the 1σ level), as can be seen from the rightmost panel. This is qualitatively similar to what we found for the globular clusters. However, the orbital eccentricity distribution for the dSph (middle panel) is very different from that of the clusters. Fewer dwarfs have very elongated orbits as their eccentricity is typically lower than 0.6. Carina even has a median eccentricity ≲ 0.2 for Models 1 and 2.

This finding leads to two interesting preliminary conclusions. Firstly, there is a weak link at most between the globular clusters in our sample and the dSph (although this is partly driven by our selection of the globular cluster sample). Secondly, the eccentricity distribution of the dwarfs is inconsistent with the predictions of cosmological simulations, where satellites are expected to be on rather radial orbits (e.g. Barber et al. 2014).

Figure 22 also confirms that Draco and Ursa Minor have very similar orbital properties and hence possibly constitute a physically connected group. Bootes I appears close to these objects in all panels ofthis figure, but as we show below, the orientation of its angular momentum differs by ~ 140 deg.

In Fig. 22 we have not plotted Leo I. This is because Leo I has extreme orbital characteristics, and is unbound in roughly 20% of the realisations for Model-1¹⁰, for example.In the cases in which a bound orbit is found, the predicted median apocentres are 819, 429, and 388 for Models 1, 2, and 3, respectively. Evidently, the apocentre of Leo I is beyond the likely virial radius of the Milky Way (estimated to be smaller than 300 kpc, see Bland-Hawthorn & Gerhard 2016, and references therein). The orbit of Leo I is quite eccentric, with medians in the range 0.6 to 0.8. The radial period is greater than 5 Gyr and its estimate varies by a factor of two for the various potentials, but a robust prediction is that the last pericentric passage (at a distance of ~ 100 kpc) took place approximately 1 Gyr ago.

The orientations of the orbital planes of the dwarfs, now including the LMC and SMC¹¹, averaged over the 10 Gyr of integration, are shown in the top panel of Fig. 23. It has been suggested that the Milky Way dwarf galaxy satellites lie on a plane (Kroupa et al. 2005). We find that their orbits tend to be perpendicular to the Galactic disc (the majority cluster at an inclination of ~ 90 ± 20 deg) but span a broad range of orientations. This implies that even though the orientation of the average plane of motion may be similar, they may rotate in the opposite sense, such as Sculptor and Ursa Minor together with Draco, the LMC and SMC (in the YZ plane). When we alternatively compare Sculptor and Sagittarius, they move in planes that are nearly perpendicular to each other (and to the Galactic disc). This ordered complexity might indicate (group) infall from a preferential direction (from a cosmic web filament aligned with the z-axis), but it seems to disfavour one single event as the cause of these configurations.

For the globular clusters we plot only the distribution of the current, instantaneous orbital plane inclination, defined as cos θ = L_z∕|L|, since the angular momentum |L| is not conserved for the majority of the clusters in our sample (for the barred potential of Model-3, L_z changes as well). The bottom panel of Fig. 23 shows that the globular clusters in our sample have a much broader distribution of inclinations than the dwarf galaxies. This behaviour might again reflect a bias in the selection of the sample of globular clusters, however, for example in the sense that there are fewer of the outer clusters that might have been associated with dwarf galaxies.

Fig. 21 Orbits integrated backward in time for 2.5 Gyr for the different dwarfs shown in different colours using the potential of Model-2.

Fig. 22 Distribution of orbital parameters for the dSph. The different colours indicate computations with different potentials. The agreement is generally good, with Model-3 systematically leading to smaller pericentres as a result of its higher mass at the radii probed by the systems. Larger differences are found for the more distant objects, revealing the sensitivity of their orbits to variations in the assumed mass distributions for the Galaxy. The dSph eccentricity distribution differs from that shown in Fig. 20 for the globular clusters in our sample. The symbols with solid error bars correspond to the median and uncertainties derived using the Monte Carlo realisations, and those with dotted error bars also take a systematic error of 0.035 mas yr−1 on each of the PM components into account.

7 Discussion

7.1 Brief summary of the systematics

The analysis performed in this paper has served to highlight the excellent quality of the Gaia DR2, and also to pinpoint its limitations. On this second aspect, we have in particular confirmed a systematic offset in the parallax of sources (see Arenou et al. 2018; Lindegren et al. 2018), which we find has an amplitude of ~ − 0.049 mas when averaged over all dSph, while averaged over the globular clusters, it is ~ − 0.025 mas. This offset varies in amplitude between the locations on the sky, which explains the difference found between our samples of dSph and globular clusters. No such offset has been found in the PMs, for example when studying open clusters (Arenou et al. 2018), from which we conclude that if such a systematic is present, it is of small amplitude. However, we do find that because of local variations (driven by the non-uniform scanning of the sky), the PMs might have an additional systematic uncertainty of ~ 0.035 mas yr⁻¹ in each direction.

7.2 Exemplifying the data quality

7.2.1 Galactic satellites

We were able to determine the PMs of all 75 globular clusters in our sample reliably (with a significance far greater than 10σ). The effect of a systematic floor noise level (see Arenou et al. 2018; Lindegren et al. 2018) is negligible for the clusters in our sample. This is also the case for the majority of the dwarf galaxies we have analysed, possibly with the exception of Leo I (because of its large distance), and Leo II (because it is small on the sky). Nonetheless, even this systematic uncertainty is typically smaller than the uncertainties of previously reported measurements of the PM of globular clusters and dSph in the Milky Way halo.

In the case of the dSph, the possibility of selecting members now via PM has revealed (and confirmed previous indications of) high spatial asymmetries, and possibly also tidal features in several systems. We did not analyse the significance of these features, nor have we attempted to establish the presence of tidal tails beyond the field of view of 2 deg that we have chosen to measure their mean PM. We expect the community to explore this area with Gaia DR2, especially now that the PMs are so accurate that reliablepredictions can be made for the orbits of these systems, and indications can be obtained on where streams, if present, might be expected.

The improved PMs for millions of stars in the Magellanic Clouds offer a unique dataset for understanding the internal dynamics of these systems. For example, Fig. 24 shows the velocity map obtained for the LMC, and reveals a high degree of order in the rotational motion of this system. This has allowed us to derive a rotation curve based on tangential velocities that is competitive to that obtained using radial velocity information. However, we also find indications in the PM residuals (after subtraction of a model of a rotating inclined disc), of streaming motion along the bar. Not only do we learn about internal structure and mass distribution of the Clouds, but this dataset will also help us in understanding how and when the two galaxies interacted, and whether and how this is related to the Magellanic stream (e.g. Kallivayalil et al. 2006a). It might even be possible to find stars stripped from the Clouds at much larger distances than previously attempted.

This brief discussion on the impressive astrometric quality of the Gaia DR2 datasets for the Galactic satellites calls for high-precision radial velocity measurements and abundances follow-up for as many of the members identified by Gaia as possible (as planned e.g. by the WEAVE and 4MOST projects, see e.g. Feltzing et al. 2017). Lists of possible members according to our analyses are given in Table D.3 for the globular clusters, dSph, and UFD galaxies, and for the Magellanic Clouds.

Fig. 23 Top: orientation of the time-averaged orbital plane, defined by the angular momentum angles (ϕ, θ) for the different dwarfs for the various potentials. The colours and error bars are the same as in Fig. 22, and open and solid circles are the median values obtained by including or excluding our estimate of the systematic error on the PMs. Most dwarfs have highly inclined orbits with respect to the Galactic plane (i.e. θ ~ 90o). Their variation in orientation (angle ϕ) over the 10 Gyr of integration is much smaller than the size of the error bars. Bottom: histogram showing the present-day inclination of the orbits for the globular clusters.

7.2.2 Substructure and debris

The methods we used thus far has relied on the objects of interest being concentrated in a specific location on the sky. Tidally torn satellites and streams, on the other hand, may extend across great parts of the celestial sphere. Here we briefly demonstrate the capacityof Gaia DR2 to also investigate this type of structure in the Milky Way using astrometric data alone. A further and deeper analysis is left to the general users of Gaia DR2.

One known substructure in the halo near the Sun that presumably originated in a disrupted satellite was discovered by Helmi et al. (1999), defining a clump in the L_z versus $\left|L_{\perp }\right|=\sqrt{L_x^2+L_y^2}$ space. The physical reality of this structure, sometimes called Helmi’s stream, was later independently confirmed by e.g. Chiba & Beers (2000); Smith et al. (2009), but no significant additional members have been identified since then. Gaia DR2 will likely reveal a very large number of new members, but measuring L_z and L_⊥ without radial velocity can only be done exactly in two small areas on the sky: in the directions of the Galactic centre and anticentre. In this case, μ_l and μ_b translate directly into the space velocities v_z and v_ϕ with knowledge of the distance D. The angular-momentum components are then simply L_z = xv_ϕ and L_⊥ = L_y ~−xv_z, where x = D + R_⊙, and R_⊙ is the Galactocentric distance of the Sun.

We show the distribution of stars within a circle of 15 degrees radius around the Galactic anticentre in the L_z vs. L_y space in Fig. 25. We considered here only stars with ϖ∕σ_ϖ∕ > 5. In addition to the dominant disc centred on (−1800, 0) kpc km s⁻¹ and the more diffuse halo centred on (0,0), stars are distinctly concentrated around (−1100, −2400) kpc km s⁻¹, corresponding to the expected location of the Helmi stream. A tight cut around the centre of this clump gives 32 candidate members (with distances from 260 pc to about 2 kpc), more than tripling the original number reported in Helmi et al. (1999). A Hertzsprung-Russel (HR) diagram of these 32 members using Gaia DR2 G magnitudes, parallaxes, and G_BP–G_RP colours is shown in Fig. 26. Despite the patchy extinction in this area of the sky, the HR diagram reveals an old, metal-poor main sequence, offset from a corresponding disc sequence by about 0.2 mag towards the blue. This figure shows four probable binaries among the 32 new members, and the three subgiants indicate a turnoff at an absolute magnitude M_G ~ 4.5. A straight extrapolation of the 32 new members from the 15-degree circle to the whole sky would give close to 2000 members in all of DR2. These can in principle be identified using radial velocity information as well.

Figure 25 shows a clear over-density only for L_y ∝ v_z < 0, and not for v_z > 0 (at fixed L_z). This stronger asymmetry in the direction of the anticentre than originally reported by Helmi et al. (1999; where the ratio was 3:1) may be used to place constraints on the accretion time and suggests that the merger may have taken place even more recently than argued by e.g. Kepley et al. (2007).

Fig. 24 PMsof stars in the LMC, represented as vectors, overlaid on a representation of the source density. This figure shows the clear consistent rotation measured by Gaia around the centre of the LMC.

Fig. 25 Distribution of DR2 stars in a circle of 15 degrees radius around the Galactic anticentre direction in the L⊥ = Ly vs. Lz space. The Helmi stream is the distinct density enhancement near the bottom of the plot (in blue). We have assumed here vrot = 220 km s−1 for the LSR velocity, vz,⊙ = 5 km s−1, and vϕ,⊙ = 7 km s−1 for the peculiar motion of the Sun, and R⊙ = 8.1 kpc.

7.3 Implications for the formation and dynamics

The unprecedentedly accurate PMs derived for the Galactic satellites from the Gaia DR2 data will allow determining the mass distribution of the Milky Way well into the realms of the dark matter halo. They will enable breaking the degeneracy between the slope of the mass density profile and the orbital anisotropy (Wilkinson & Evans 1999; Watkins et al. 2010), the latter being the limiting factor thus far, which Gaia DR2 has turned into an observable. Interestingly, our measurements indicate that the orbits of the dSph are not very radial, and this appears to challenge expectations derived from cosmological simulations in the ΛCDM framework (e.g. Lux et al. 2010; Cautun & Frenk 2017). Their relatively low eccentricity (r_apo∕r_peri < 4) and non-penetrating orbits (r_peri > 30 kpc) may also disfavour models for the transition of dwarf irregulars into dSph via a tidal-stirring mechanism (Kazantzidis et al. 2017).

We have found relatively small differences in the orbits in three different realistic Galactic potentials when we integrated over short timescales (the potentials are based on Allen & Santillan 1991; Robin et al. 2003; McMillan 2017). This is the case for globular clusters that do not probe the inner few kpc of the Milky Way and for the dwarf galaxies. This implies that the orbit-forecasting power is high, and it might therefore be used to search for tidal tails particularly for the globular clusters, both through their predicted location on the sky but also using PM information that is now available thanks to Gaia DR2. This is a particularly interesting avenue because of the very high constrainingpower of streams on the mass distribution in our Galaxy (see Johnston & Carlberg 2016, and references therein).

Anotherinteresting constraint comes from the most distant satellite Leo I. This object has most recently been used to derive a limit on the mass of the dark matter halo of the Milky Way on the basis of the so-called timing argument by Sohn et al. (2013). From our own measurements of the PM of Leo I, we find that it is barely bound, with its orbit extending well beyond recent estimates of the virial radius of our Galaxy. Although it is very unlikely that Leo I is unbound (Boylan-Kolchin et al. 2013), we may use it to derive a lower limit to the mass of the Milky Way if we assume it has the escape velocity. Neglecting the contribution of the disc(s) and bulge, for a (non-truncated) NFW halo, $$M_{MW}(r_{\text{LeoI}})=\frac{v_{\text{LeoI}}^2{r}_{\text{LeoI}}}{2G}\frac{\log (1+x_{\text{LeoI}})-x_{\text{LeoI}}/(1+x_{\text{LeoI}})}{\log (1+x_{\text{LeoI}})},$$(6)

with x_LeoI = r_LeoI∕r_s. When we use our estimates of $v_{\text{LeoI}}~{217.3}_{-48.9}^{+62.6}$ km s⁻¹ and $r_{\text{LeoI}}~{257.8}_{-35.1}^{+16.8}$ kpc, this implies a lower limit for the enclosed mass of the Milky Way of $M_{MW}(r_{\text{LeoI}})={9.1}_{-2.6}^{+6.2}\times {10}^{11}{M}_{\odot }$ assuming r_s = 18.6 kpc as in McMillan (2017), where the error bars indicate the 16th and 84th percentiles. When we let this parameter vary in the range 10 ≤ r_s ≤ 30 kpc, the estimate of the lower limit to the virial mass varies by ~ 15%, that is, within the uncertainties bracketed by the measurement errors. This value is in line with previous work (see review by Bland-Hawthorn & Gerhard 2016, and references therein) and does not preclude a light dark matter halo for our Galaxy (of ~ 10¹² M_⊙, e.g. Battaglia et al. 2005; Gibbons et al. 2014).

Leo I is an intriguing object because it is very distant and its velocity indicates that it is receding from us. A possible explanation for a system such as Leo I is that it has experienced a three-body interaction with the Magellanic Clouds (Sales et al. 2007). Although we do not explore this possibility here, we find that the phase-space distribution of dwarf galaxies is not homogeneous. These satellites tend to have orbits with angular momenta perpendicular to the Galactic disc, meaning that their orbits take place in planes with varying orientations, but always with high inclination. This supports the idea of filamentary infall (see Libeskind et al. 2005), and might also imply that some of the dSph in our sample have fallen in together as a group (Lynden-Bell & Lynden-Bell 1995; Li & Helmi 2008). For example, the relative distances of Ursa Minor and Draco remain relatively small when computed over our 10 Gyr long orbital integrations, thus favouring some amount of group infall.

Our measurements rule out that the dwarf galaxies are on one single narrow “disc”, an idea that has led to an important debate in the literature. Furthermore, although some objects are on the same plane, they rotate in a different sense around the Galaxy, making it less likely that they formed in single event such as a major merger (Yang et al. 2014, this scenario would also require that their eccentricities be relatively high, in contrast to what we find from Gaia DR2 data). On the other hand, our measurements will finally enable exploring how the satellite population was put in place and how this relates to the environment of the Milky Way on a firm basis (Libeskind et al. 2015).

Future studies of the orbital properties of the more distant globular clusters will allow us to establish their relation to the present-day dwarf galaxy population (and further test the ideas put forward in Lynden-Bell & Lynden-Bell 1995). On the other hand, some of the globular clusters in our sample may have been associated with long-gone accreted galaxies, whose debris we might expect will be discovered and characterised using Gaia DR2 data. Together, this will shed light on the build-up of the globular cluster population and its link to former and current satellite galaxies.

Fig. 26 Hertzprung-Russell diagram using Gaia DR2 G magnitudes, colours, and parallaxes in the direction of the Galactic anticentre. The yellowish cloud represents a random sample of disc stars, selected on the basis of their disc-like Lz and L⊥, while the blue dots are the 32 new members of the Helmi stream.

8 Conclusions

The second data release from the Gaia mission has delivered its promise of new and accurate PM measurements for a billion stars across the full Galaxy and its nearest neighbours. Our analysis of the PMs of stars in roughly half of the Galactic globular clusters, in all the known dSph galaxies, and in the Magellanic Clouds has allowed us to derive their mean motions as they orbit the Galaxy much more precisely than ever before, despite systematics that are clearly present in this data release. The simple analyses carried out in this paper have confirmed, and also revealed for many globular clusters, previously reported internal dynamical complexity (such as rotation and the presence of extended halos). The PMs, and hence the orbits, of the dSph have finally been pinned-down, and this has uncovered their relatively coherent phase-space distribution, which is not consistent with a single “disc of satellites”, however. The astounding dynamical maps of the Magellanic Clouds contain such richness that is has even been possible to derive for the first time a high quality rotationcurve for the Large Magellanic Cloud based on tangential velocities alone and to unveil the dynamical imprint of the bar.

Much remains to be understood and discovered from the PMs that have become available with Gaia DR2. This dataset will undoubtedly keep the Galactic astronomy community busy for many years to come.

Acknowledgements

This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia. The Gaia mission and data processing have financially been supported, in alphabetical order by country, by the Algerian Centre de Recherche en Astronomie, Astrophysique et Géophysique of Bouzareah Observatory; the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (FWF) Hertha Firnberg Programme through grants T359, P20046, and P23737; the BELgian federal Science Policy Office (BELSPO) through various PROgramme de Développement d’Expériences scientifiques (PRODEX) grants and the Polish Academy of Sciences - Fonds Wetenschappelijk Onderzoek through grant VS.091.16N; the Brazil-France exchange programmes Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Coordenação de Aperfeicoamento de Pessoal de Nível Superior (CAPES)-Comité Français d’Évaluation de la Coopération Universitaire et Scientifique avec le Brésil (COFECUB); the Chilean Dirección de Gestión de la Investigación (DGI) at the University of Antofagasta and the Comité Mixto ESO-Chile; the National Science Foundation of China (NSFC) through grants 11573054 and 11703065; the Czech Republic Ministry of Education, Youth, and Sports through grant LG 15010, the Czech Space Office through ESA PECS contract 98058, and Charles University Prague through grant PRIMUS/SCI/17; the Danish Ministry of Science; the Estonian Ministry of Education and Research through grant IUT40-1; the European Commission’s Sixth Framework Programme through the European Leadership in Space Astrometry (https://www.cosmos.esa.int/web/gaia/elsa-rtn-programme, ELSA) Marie Curie Research Training Network (MRTN-CT-2006-033481), through Marie Curie project PIOF-GA-2009-255267 (Space AsteroSeismology & RR Lyrae stars, SAS-RRL), and through a Marie Curie Transfer-of-Knowledge (ToK) fellowship (MTKD-CT-2004-014188); the European Commission’s Seventh Framework Programme through grant FP7-606740 (FP7-SPACE-2013-1) for the Gaia European Network for Improved data User Services (https://gaia.ub.edu/Twiki/bin/view/GENIUS/WebHome, GENIUS) and through grant 264895 for the Gaia Research for European Astronomy Training (https://www.cosmos.esa.int/web/gaia/great-programme, GREAT-ITN) network; the European Research Council (ERC) through grants 320360 and 647208 and through the European Union’s Horizon 2020 research and innovation programme through grants 670519 (Mixing and Angular Momentum tranSport of massIvE stars – MAMSIE) and 687378 (Small Bodies: Near and Far); the European Science Foundation (ESF), in the framework of the Gaia Research for European Astronomy Training Research Network Programme (https://www.cosmos.esa.int/web/gaia/great-programme, GREAT-ESF); the European Space Agency (ESA) in the framework of the Gaia project, through the Plan for European Cooperating States (PECS) programme through grants for Slovenia, through contracts C98090 and 4000106398/12/NL/KML for Hungary, and through contract 4000115263/15/NL/IB for Germany; the European Union (EU) through a European Regional Development Fund (ERDF) for Galicia, Spain; the Academy of Finland and the Magnus Ehrnrooth Foundation; the French Centre National de la Recherche Scientifique (CNRS) through action “Défi MASTODONS”, the Centre National d’Études Spatiales (CNES), the L’Agence Nationale de la Recherche (ANR) “Investissements d’avenir” Initiatives D’EXcellence (IDEX) programme Paris Sciences et Lettres (PSL*) through grant ANR-10-IDEX-0001-02, the ANR “Défi de tous les savoirs” (DS10) programme through grant ANR-15-CE31-0007 for project “Modelling the Milky Way in the Gaia era” (MOD4Gaia), the Région Aquitaine, the Université de Bordeaux, and the Utinam Institute of the Université de Franche-Comté, supported by the Région de Franche-Comté and the Institut des Sciences de l’Univers (INSU); the German Aerospace Agency (Deutsches Zentrum für Luft- und Raumfahrt e.V., DLR) through grants 50QG0501, 50QG0601, 50QG0602, 50QG0701, 50QG0901, 50QG1001, 50QG1101, 50QG1401, 50QG1402, 50QG1403, and 50QG1404 and the Centre for Information Services and High Performance Computing (ZIH) at the Technische Universität (TU) Dresden for generous allocations of computer time; the Hungarian Academy of Sciences through the Lendület Programme LP2014-17 and the János Bolyai Research Scholarship (L. Molnár and E. Plachy) and the Hungarian National Research, Development, and Innovation Office through grants NKFIH K-115709, PD-116175, and PD-121203; the Science Foundation Ireland (SFI) through a Royal Society - SFI University Research Fellowship (M. Fraser); the Israel Science Foundation (ISF) through grant 848/16; the Agenzia Spaziale Italiana (ASI) through contracts I/037/08/0, I/058/10/0, 2014-025-R.0, and 2014-025-R.1.2015 to the Italian Istituto Nazionale di Astrofisica (INAF), contract 2014-049-R.0/1/2 to INAF dedicated to the Space Science Data Centre (SSDC, formerly known as the ASI Sciece Data Centre, ASDC), and contracts I/008/10/0, 2013/030/I.0, 2013-030-I.0.1-2015, and 2016-17-I.0 to the Aerospace Logistics Technology Engineering Company (ALTEC S.p.A.), and INAF; the Netherlands Organisation for Scientific Research (NWO) through grant NWO-M-614.061.414 and through a VICI grant (A. Helmi) and the Netherlands Research School for Astronomy (NOVA); the Polish National Science Centre through HARMONIA grant 2015/18/M/ST9/00544 and ETIUDA grants 2016/20/S/ST9/00162 and 2016/20/T/ST9/00170; the Portugese Fundação para a Ciência e a Tecnologia (FCT) through grant SFRH/BPD/74697/2010; the Strategic Programmes UID/FIS/00099/2013 for CENTRA and UID/EEA/00066/2013 for UNINOVA; the Slovenian Research Agency through grant P1-0188; the Spanish Ministry of Economy (MINECO/FEDER, UE) through grants ESP2014-55996-C2-1-R, ESP2014-55996-C2-2-R, ESP2016-80079-C2-1-R, and ESP2016-80079-C2-2-R, the Spanish Ministerio de Economía, Industria y Competitividad through grant AyA2014-55216, the Spanish Ministerio de Educación, Cultura y Deporte (MECD) through grant FPU16/03827, the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia “María de Maeztu”) through grant MDM-2014-0369, the Xunta de Galicia and the Centros Singulares de Investigación de Galicia for the period 2016-2019 through the Centro de Investigación en Tecnologías de la Información y las Comunicaciones (CITIC), the Red Española de Supercomputación (RES) computer resources at MareNostrum, and the Barcelona Supercomputing Centre - Centro Nacional de Supercomputación (BSC-CNS) through activities AECT-2016-1-0006, AECT-2016-2-0013, AECT-2016-3-0011, and AECT-2017-1-0020; the Swedish National Space Board (SNSB/Rymdstyrelsen); the Swiss State Secretariat for Education, Research, and Innovation through the ESA PRODEX programme, the Mesures d’Accompagnement, the Swiss Activités Nationales Complémentaires, and the Swiss National Science Foundation; the United Kingdom Rutherford Appleton Laboratory, the United Kingdom Science and Technology Facilities Council (STFC) through grant ST/L006553/1, the United Kingdom Space Agency (UKSA) through grant ST/N000641/1 and ST/N001117/1, as well as a Particle Physics and Astronomy Research Council Grant PP/C503703/1. We are verygrateful to the referee, Mario Mateo, for his prompt, very positive, and constructive report. We are grateful to Dana Casetti-Dinescu, Mike Boylan-Kolchin and Gary Mamon for suggesting various small corrections.

Appendix A Additional descriptions

A.1 Example of data extraction

To illustratethe different stages of the data extraction, we show here an example following the procedure used for globular cluster NGC 2298.

Fig. A.1 PM of stars within 0.5 deg. from the centre of NGC2298 as extracted from the GACS archive. Left: all ~ 259 000 stars. Middle: zoom-in on the PM field, which shows a concentration of the stars near position (3.2, −2.2) mas yr−1 corresponding to cluster members. Right: another zoom-in on the PM field for the cluster.

Fig. A.2 From left to right: distribution of positions on the sky for the first selection on cluster star PMs, showing the cluster in the centre, PM distribution for the cluster field, and distribution of parallax vs. PM in RA.

A first selection was made in a field with a radius of about 1 degree around the cluster centre as given in Harris10. The PMs for ~ 259 000 stars in that field are shown in the left panel of Fig. A.1. A zoom-in on the diagram (middle panel) shows a small concentration in the PM field near position (3.2, −2.2) mas yr⁻¹, which is likely to be the cluster. No other such concentrations are apparent in the same field. A further zoom-in on the concentration (right panel of Fig. A.1) shows it to be real. The data in this field are selected as a new subset, for which the distribution of positions on the sky are shown in the left panel of Fig. A.2, which reveals the cluster more clearly. The PM distribution for this new coverage is also improved, showing much reduced field star contamination (middle panel). The distribution over parallaxes now also clearly shows the cluster with little contamination (right panel). Finally, for the pre-selection procedure, the colour-magnitude diagram was created for the selection (Fig. A.3), which now contains about 1500 stars. The selected data thus obtained were saved as a CSV file for further analysis through the cluster analysis software.

The analysis software iterates over an astrometric solution for the cluster, in which the parallax and PM of the cluster are solved for. The cluster centre is determined based on the mean position of all selected members, without any weighting being applied. First the cluster stars positions are projected on a tangential plane with the assumed cluster centre as zero point. The mean of the Cartesian positions is determined, which then is de-projected to provide the new cluster centre on the sky. The relevant equations have been presented in Gaia Collaboration (2017).

Fig. A.3 HR diagram for the cluster field before (left) and after (right) the astrometric solution.

The astrometric solution is based on the parallax and PM determinations for the individual member stars. For each star the contributionto the solution was normalised by the covariance matrix (for more details, see Appendix A.1 in Gaia Collaboration 2017). This created for each star three uncorrelated contributions with unit-weight error variance, which contribute to the cluster astrometric parameter solution. Outliers were rejected on the basis of the normalised residuals, except in a few cases where field stars significantly disturbed the astrometric parameters. In these cases, which were mainly found at lower galactic latitudes, a 3σ cut was applied to the individual parallaxes and PMs as based on their standard uncertainties and the estimate of the cluster PM and parallax (see also Fig. A.4). The astrometric solution removes a large fraction, but not all, of the remaining field stars, as canbe seen from the HR diagram shown in the right panel of Fig. A.3. When this residual population of field stars still significantly disturbed the HR diagram, we performed a final cleanup in TOPCAT on the basis of that diagram.

Fig. A.4 Two examples of the parallax and PM dispersion diagrams for clusters at a distance of approximately 2 (top) and 5 (bottom) kpc. The parallax and PMs are plotted against the standard uncertainties. The blue lines show the expected uncertainty levels relative to the mean value at –2, –1 0 1 and 2 σ, including noise contributions from the internal velocity dispersion and the dispersion of the parallaxes from the depth of the cluster. At low levels of standard uncertainty, the contributions from the internal dispersion are easily detected. The colouring ofthe data points reflects the error correlations: red means strongly negative, blue means strongly positive, and black means low.

Two criteria are set in the iterations over the astrometric solution, one of which could require an extended extract from the archive. The first criterion concerns the internal PM dispersion. This could be detected from the dispersion in the PMs of the brightest cluster members. This is, however, not an unambiguous process, as the observed internal velocity dispersion is often observed to be strongly dependent on distance from the cluster centre. The dispersions were examined in the radial and transverse directions, and three examples are shown in Fig. A.5. The modelling of these dependencies is, however, beyond the scope of the present paper, and an average dispersion was used in all cases for the whole cluster. The values found were generally between 2 and 8 km s⁻¹.

Strong variations in coverage of stars with five-parameter solutions in the astrometry (Lindegren et al. 2018) also show the relation between standard uncertainties on the astrometric parameters and the brightness of the stars. With relatively homogeneous coverage, the distribution is narrow, while, as is the case for ω Cen, a multi-layered set of relations is observed for poorly covered objects (see Fig. A.6). An additional feature appears to be the increased noise on the astrometric data for variable stars with large amplitudes, such as the RR Lyrae stars in globular clusters, as shown in Fig A.7. The exact reason for this is still unclear.

Fig. A.5 Three examples of the PM dispersion in the radial and transverse directions as a function of distance to the cluster centre. The three vertical lines indicate, from left to right, the core radius, the half-light radius, and the tidal radius following Harris10.

Fig. A.6 Two examples of astrometric data coverage with five-parameter solutions. Left: ω Cen, the worst case; right: NGC 5272, a more average example of coverage. The gaps in the coverage for ω Cen are the result of the filters that have been applied to the astrometric data. The cyan circles are at intervals of 35 pc in ω Cen and 10 pc in NGC 5272.

A.2 Exemplifying some systematics in Gaia DR2 using dwarf galaxies

Figure A.8 shows the parallax zero-point offset for the stars in Carina and demonstrates that it is very reliably measured and is independent of the magnitude of the stars used. This indicates a systematic effect present in the parallax measurements of the DR2 data, which for Carina is of the order of ϖ_DR2 − ϖ_lit = −0.015 − 0.0095 ~ − 0.024 mas.

Figure A.9 shows the gridding pattern present in the parallaxes, in this case, for stars in the LMC and the SMC (see also Lindegren et al. 2018). This pattern has an amplitude of ~ 0.03 mas (see the bottom panel of Fig. A.10), and is clearly apparent when analysing sufficiently large objects on the sky, but it is likely to be present throughout the full sky (see e.g. Arenou et al. 2018).

Fig. A.7 From left to right: standard uncertainties on the astrometric parameters of individual stars as a function of magnitude in ω Cen and NGC 5272, and standard uncertainties on the photometric data for NGC 5272. The photometric data clearly show the variability of the RR Lyrae stars (around G = 15.3), while the astrometric data also show locally much increased uncertainties at the same brightness.

Fig. A.8 Parallax of the stars in the Carina dSph as function of their G magnitude. A clear systematic offset is apparent and is significantly measured. The expected parallax (based on literature values) is indicated by the red dashed line.

Figure A.10 shows the variation in PM and parallax along stripes of width 0.2° in the LMC after a model for the PM field has been subtracted. The variation seen is produced both by the systematic errors in Gaia DR2 and the shortcomings of the simple model.

The correlation between the PMs of individual stars in the field of view towards the dSph in our sample is is illustrated in Fig. A.11, where we have plotted the individual stars with different colours thatindicate the amount of correlation (from very negative, to none, to very positive) between ${\mu }_{\alpha }^{*}$ and μ_δ. Fornax in particular shows a strong correlation in the PM components of the individual stars, and these are of course then reflected in the PM correlation coefficient given in the table. The amplitude and orientation of the correlation differs among the dwarf galaxies, indicating that the correlations are localised on the sky and do not have the same amplitude everywhere. Furthermore, there are regions where these correlations are negligible when averaged out, as in the case of Carina, Draco, Sculptor, and even the Sagittarius dSph. We also note that all the stars in the field of view towards these dwarfs are affected in a similar way, and not only the members.

Fig. A.9 Measured parallaxes of stars in the LMC (upper panel) and SMC (lower panel). The banding associated with the Gaia scanning law is clearly visible.

Fig. A.10 Parallaxes and residual PMs of sources in the LMC. The lines show the median value as a function of coordinate position x for all LMC sources within 0.2° stripes centred at y = −0.5, 0 or 0.5°. Residual PMs are calculated after subtraction of a model disc with parameters determined using only sources within 3° of the centre. The non-zero values are due to both the systematic errors in the Gaia data and the differences between the simple disc model and the true dynamics of the LMC.

Fig. A.11 Correlations in the PMs of the stars in the field of view towards the different dSph galaxies in our sample. The different colours indicate the amplitude of the correlations.

Appendix B Details of the LMC and SMC modelling

In this appendix we provide a rigorous description of the coordinate system used in Sects. 2.2 and 5, the modelling assumptions used to derive Eq. (3), and the deprojection of the PMs shown in Fig. 18. We refer to the centre of the Cloud as C, to the observer as O, and to a source in the Cloudas S.

The vectors [x y z] form an inertially fixed right-handed orthogonal triad with z the unit vector from observer to C at the reference epoch, x the unit vector in the direction of increasing α at C, and y the unit vector in the direction of increasing δ at C.

B.1 Position and PM in the xyz system

In terms of the celestial position (α, δ) and PM components (μ_α*, μ_δ), we have the unit vector from the observer to a source, u, given by $$u=r,\dot{u}=p{\mu }_{\alpha *}+q{\mu }_{\delta },$$(B.1)

where [p q r] is the local normal triad at (α, δ) $$p=\left[\begin{array}{c}-\sin \alpha \\ \cos \alpha \\ 0\end{array}\right],q=\left[\begin{array}{c}-\sin \delta \cos \alpha \\ -\sin \delta \sin \alpha \\ \cos \delta \end{array}\right],r=\left[\begin{array}{c}\cos \delta \cos \alpha \\ \cos \delta \sin \alpha \\ \sin \delta \end{array}\right].$$(B.2)

The xyz system coincides with the local normal triad at (α_C, δ_C) $$x=\left[\begin{array}{c}-\sin {\alpha }_C\\ \cos {\alpha }_C\\ 0\end{array}\right],y=\left[\begin{array}{c}-\sin {\delta }_C\cos {\alpha }_C\\ -\sin {\delta }_C\sin {\alpha }_C\\ \cos {\delta }_C\end{array}\right],z=\left[\begin{array}{c}\cos {\delta }_C\cos {\alpha }_C\\ \cos {\delta }_C\sin {\alpha }_C\\ \sin {\delta }_C\end{array}\right].$$(B.3)

The components of $u,\dot{u}$ in the xyz system are obtained as scalar products (e.g. x = x×u, $\dot{x}=x\times \dot{u}$) from which we can derive Eqs. (2), where we refer to ẋ, ẏ as μ_x, μ_y (they are not strictly speaking PMs, but it is convenient to give them this notation).

When (α_C, δ_C) are used as a fixed reference point, the Cartesian coordinates (x, y, z, ẋ, ẏ, ż) provide a useful substitute for (α, δ, μ_α*, μ_δ). The six components are a redundant set, and when working in a limited area around C (in principle as long as z> 0, i.e. within 90° from C), it is possible to use the non-redundant set (x, y, ẋ, ẏ), with $$z=\sqrt{1-x^2-y^2},\dot{z}=-(x\dot{x}+y\dot{y})/z.$$(B.4)

(x, y) is equivalent to the orthographic projection in cartography.¹²

B.2 Kinematic model

Assuming a flat disc, the vector R from C to S must be in the plane of the disc, which gives the condition $$n\times R=0,$$(B.5)

where we have reintroduced the vector n, the normal to the disc plane (such that rotation about n is positive), and we also reintroduce the two normal unit vectors l and m, which form a right-handed triad with n (e.g. Eq. (4)). A source in the plane of the disc can be described in terms of rectangular coordinates ξ, η, where $$R=l\xi +m\eta .$$(B.6)

The motion of S is the vectorial sum of the bulk motion (of C) and the peculiar motion of S with respect to C. If we assume that the peculiar motion is circular with angular velocity ω(R), we have $$v_S=v_C+(m\xi -l\eta )\omega (R),$$(B.7)

where v_C is the bulk motion and $R=\sqrt{{\xi }^2+{\eta }^2}$.

At the reference epoch, we can write R in terms of the position of the source s and the centre c as R = s−c = us −z, where s is the distance to the source. Introducing the inverse distance factor f = s⁻¹ and inserting in Eq. (B.5) gives $$f=(n\times u)/n_z=ax+by+z={(1-x^2-y^2)}^{1/2}+ax+by,$$(B.8)

where $$a=n_x/n_z,b=n_y/n_z,$$(B.9)

and n_x = x×n etc are the components of n in the xyz system. These components, along with those of l and m and the values a and b are given as a function of the inclination i and the line-of-nodes position angle Ω in the main text (Eqs. (4) and (5)). We now find $$\begin{array}{ll}\xi \hfill & =l\times R=\frac{l_{x}x+l_{y}x}{z+ax+by}\hfill \\ \eta \hfill & =m\times R=\frac{(m_x-a{m}_z)x+(m_y-b{m}_z)y}{z+ax+by}.\hfill \end{array}$$(B.10)

Turning now to the PMs, we seek the corresponding relationships between (ẋ, ẏ) and $(\dot{\xi },\dot{\eta })$. We know that $$\dot{\xi }=-\eta \omega (R),\dot{\eta }=\xi \omega (R).$$(B.11)

The PM vector is $$\dot{u}=\frac{\text{d}(sf)}{\text{d}t}=\left(v_S-u(u\times v_S)\right)f,$$(B.12)

which, using our previous results, we can rewrite as $$\begin{array}{l}\dot{u}=\left(v_C-u(u\times v_C)\right)(ax+by+z)+((m-u(u\times m))(l_{x}x+l_{y}y)-\\ (l-u(u\times l))((m_x-a{m}_z)x+(m_y-b{m}_z)y))\omega (R).\end{array}$$(B.13)

Taking the scalar products with x and y gives explicit expressions for ẋ and ẏ as functions of x and y. If ω is constant, the expressions contain terms up to the third power in x and y. At C (x = y =0) we find, as expected, $$\dot{x}=v_x,\dot{y}=v_y,$$(B.14)

where v_x, v_y, v_z, etc. are the components of v_C in the xyz system.

Retaining only first-order terms in x and y while assuming constant ω, and then taking derivatives, we have $$\begin{array}{ll}\partial \dot{x}/\partial x\hfill & =a{v}_x-v_z+a{l}_x{m}_z\omega ,\hfill \\ \partial \dot{x}/\partial y\hfill & =b{v}_x-n_z\omega +b{l}_x{m}_z\omega ,\hfill \\ \partial \dot{y}/\partial x\hfill & =a{v}_y+n_z\omega +a{l}_y{m}_z\omega ,\hfill \\ \partial \dot{y}/\partial y\hfill & =b{v}_y-v_z+b{l}_y{m}_z\omega ,\hfill \end{array}$$(B.15)

which hold exactly at C if d ω(r)∕dR = 0 for R = 0. These equations (writing ẋ, ẏ as μ_x, μ_y) were used in this study to determine v_z, i, Ω, and ω.

The orientation of the LMC plane is given by the unit vector n, the direction of which is conventionally given by the two angles i and Ω: $$n=x\sin i\cos \Omega -y\sin i\sin \Omega +z\cos i,$$(B.16)

from which (with the definition of l and m) we have the components of all three of these vectors in the xyz system, as given in Eq. (4).

B.3 Estimating the kinematic parameters for fixed C and constant ω

From Eq. (B.15) we see that (ẋ, ẏ) should vary linearly with (x, y) for constant ω. This is a reasonable approximation in the inner few degrees of the LMC. Fitting the linear relation $$\begin{array}{ll}\dot{x}\hfill & =v_x+A_{x}x+A_{y}y\hfill \\ \dot{y}\hfill & =v_y+B_{x}x+B_{y}y\hfill \end{array}$$(B.17)

we immediately obtain estimates of v_x, v_y, and the four gradients ∂ẋ∕∂x = A_x, etc. We can use Eq. (B.15) to express these gradients as functions of the six kinematic parameters v_x, v_y, v_z, i, Ω, and ω (assumed constant). In this study we usually then held v_z constant, in which case we varied i, Ω, and ω to minimise the sum of the square residuals $$S={(\partial \dot{x}/\partial x-A_x)}^2+{(\partial \dot{x}/\partial y-A_y)}^2+{(\partial \dot{y}/\partial x-B_x)}^2+{(\partial \dot{y}/\partial y-B_y)}^2,$$(B.18)

where A_x etc. are measured from the data, and ∂ẋ∕∂x are predictions of the model. We can also leave v_z free, in which case Eqs. (B.15) can be directly solved to determine v_z, i, Ω, and ω.

B.4 De-projection method

The observed PMs, or ẋ, ẏ, are the projections of the true space motions on the celestial sphere, or normal to z. This projection from 3D to 2D cannot be inverted, but if we assume that the true motions are confined to the (known) plane of the LMC, itis possible to project the observed motions back to that plane. We call this de-projection.

We assume that v_C, i, and Ω are known. The starting point is Eq. (B.12), but with the velocity of S written as $$v_S=v_C+l\dot{\xi }+m\dot{\eta },$$(B.19)

which allows arbitrary motions in the plane. We have then $$\dot{u}=\left(v_C-u(u\times v_C)\right)f+(l-u(u\times l))f\dot{\xi }+(m-u(u\times m))f\dot{\eta }.$$(B.20)

Taking the scalar products with x and y gives two linear equations, $$\begin{array}{ll}(l_x-x(l_{x}x+l_{y}y))\dot{\xi }\hfill & +(m_x-x(m_{x}x+m_{y}y+m_{z}z))\dot{\eta }\hfill \\ \hfill & =(\dot{x}-v_x+x(v_{x}x+v_{y}y+v_{z}z))/(ax+by+z)\hfill \\ (l_y-y(l_{x}x+l_{y}y))\dot{\xi }\hfill & +(m_y-y(m_{x}x+m_{y}y+m_{z}z))\dot{\eta }\hfill \\ \hfill & =(\dot{y}-v_y+y(v_{x}x+v_{y}y+v_{z}z))/(ax+by+z)\hfill \end{array}$$(B.21)

from which $\dot{\xi }$ and $\dot{\eta }$ can be solved. The corresponding position (ξ, η) is obtained from (B.10).

To map the kinematics of the LMC, it is more convenient to transform the Cartesian $(\xi ,\eta ,\dot{\xi },\dot{\eta })$ into polar coordinates R, ϕ and the corresponding velocity components v_R (in the direction of increasing R) and v_T (tangential velocity), as in Fig. 18.

Appendix C Globular clusters and dSph solution data

Appendix D Orbital integrations

Figure D.1 compares the circular velocity curves for the three different Galactic potentials. The curves are quite similar between ~ 3 and ~ 40 kpc, but they differ substantially in the inner as well as in the outer Galaxy. This will lead to some of the differences in the orbits, as discussed in Sect. 6.

Figure D.1 shows some of these differences for the subset of globular clusters shown in Fig. 19. We plot here the orbits in cylindrical coordinates and for a shorter period of time to show more clearly how the orbits diverge from each other due to the differences in the gravitational potentials of the various models.

Figure D.3 shows the distribution of orbital parameters for the dSph in our sample, derived by drawing 1000 Monte Carlo realisations of the observables taking into account their uncertainties. The correlations seen are not all due to correlations in the errors themselves, but also reflect that orbital parameters are not really fully independent. It is important to bear this in mind when interpreting Figs. 20 and 22, where we have plotted uncorrelated error bars to facilitate visual inspection.

Fig. D.1 Circular velocity curves for the three Galactic potentials considered for the orbit integrations in Sect. 6.

Fig. D.2 Orbits of the globular clusters shown in Fig. 19, now plotted in cylindrical coordinates and as a function of time (right panels) for shorter integration times (approximately three radial oscillations). The different colours correspond to the axisymmetric potentials of Models 1 (blue) and 2 (black) and to the barred Model-3 (red). The solid circles in the leftmost panels denote the present-day positions, while the squares are the positions for each of the models at the end of the chosen integration time.

Fig. D.3 Distribution of orbital parameters for the dwarfs (in different colours) showing the results obtained from the different Monte Carlo realisations (only for those within the 1σ uncertainties on the observables), computed using the potential of Model-1.

References

All Tables

All Figures

Fig. 1 Top: distance distribution of the 75 globular clusters included in the present study. Bottom: standard uncertainties on the PM in declination as a function of number of cluster members nMemb used in the solution. The diagonal line represents a fit to the relation $\sigma {\mu }_{\delta }=a/\sqrt{nMemb}$ where we find a = 0.3 [mas/yr]. A similar dependence on the number of members is found for the parallax uncertainty (with a = 0.15 [mas]) and for σμα* (where a= 0.25 [mas/yr]).

In the text

Fig. 2 Field-of-view towards the dSph galaxies (the nine classical and one UFD) in our sample. The stars shown correspond to membersaccording to the photometric selection (on the RGB and BHB) and the astrometric procedure (within 2σ from the mean PM of the object). The striping apparent in Sagittarius and Sextans is driven in part by the scanning law. The hole in the centre of Sagittarius corresponds to the location of the globular cluster NGC 6715 (M54).

In the text

Fig. 3 Colour-magnitude of the stars in the field of view towards the dSph galaxies (the eight classical and one UFD) in our sample. The blue lines mark our (relatively tight) pre-selection of tentative members (on the RGB and BHB) that is fed to the pipeline to derive mean PMs. The coloured points indicate stars within 3σ of our determination of the mean PM of the object. This means that cyan points satisfy both the PM and CMD selections.

In the text

	Fig. 4 Comparison between the tidal radii rt (according to Harris10) of the 75 globular clusters in our sample and the maximum radii at which we have been able to detect cluster members in the present study. The diagonal line represents the one-to-one relation.
In the text

	Fig. 5 Density distribution on the sky of the stars selected as members of the LMC (left) and SMC (right). Positions are shown in the (x, y) coordinatesdescribed in the text (Eq. (2)). In all figures that use this coordinate system, the x-axis has been inverted so that it corresponds to the usual inversion of right ascension.
In the text

	Fig. 6 Comparison between parallaxes as derived from the Gaia DR2 data and parallaxes derived from the cluster distances as given in Harris10.
In the text

Fig. 7 Top: parallax error against the Gaia parallax as determined from the Gaia data for 75 globular clusters. The black curves are the 1 and 3σ limits. The vertical dashed line corresponds to 3 times our estimate of the systematic error on the parallax. Bottom: PM errors against the PMs in right ascension (open circles) and declination (solid circles) for the clusters in our sample. The curves represent the value of the PM for 100σμ. The PM measurement has a significance lower than 10σ only for NGC 6453, for which ${\mu }_{\alpha *}/{\sigma }_{{\mu }_{\alpha *}}~4$.

In the text

	Fig. 8 Comparison of the Gaia PMs (in right ascension: top, and declination: bottom) to measurements reported in Dinescu et al. (1999, 2003); Casetti-Dinescu et al. (2007, 2010, 2013) for 31 globular clusters.
In the text

	Fig. 9 Top: comparison between the ground-based and Gaia radial velocities for 52 clusters with at least two measurements. Bottom: standard uncertainties on the mean radial velocities as a function of the number of stars contributing to the mean. The diagonal line shows the effect of an additional contribution of 4 km s−1 originating from the internal velocity dispersion.
In the text

	Fig. 10 Sky distribution of the 75 globular clusters in our sample in Galactic coordinates. Their tangential velocities are denoted by the size and direction of the arrows, while the colours indicate their line-of-sight velocities. The inset shows a zoom-in of the central 60 × 60 deg2.
In the text

Fig. 11 PM “systematics” in radial (blue triangles) and transverse (red squares) directions, as a function of distance to the cluster centre. The cyan line is a fit to the variation of the mean radial component of the PMs as a function of distance, while the grey line is the expected trend resulting from perspective contraction or expansion. From left to right: NGC 6656 shows a strong rotation signal as well as perspective expansion; NGC 3201 shows no rotation but very strong perspective contraction; and NGC 6397 shows halo expansion associated with core collapse, clearly different from the signal expected from the small perspective contraction. The vertical line indicates the tidal radius of the cluster reported in Harris10.

In the text

Fig. 12 PMs of the stars in the field of view towards the different dSph galaxies in our sample. Members, defined as stars within 3σ of the mean measured PM and located in the expected region of the CMD, are shown in cyan. The green points correspond to stars that also fall in the CMD-selected box, but are not within 3σ of the systemic PM. Especially for Sculptor and Fornax, it is quite clear that there may be more members, but very likely, the large errors on the PMs of individual stars place them beyond the 3σ ellipse.

In the text

	Fig. 13 Counts in 0.2 deg wide bins on the sky with at least 100 stars for stars in the Sagittarius dwarf (left) and the average PM in μα* (middle) and μδ (right) for each of these bins.
In the text

	Fig. 14 Blue histograms: distribution of average PM in the bins shown in Fig. 13; red histograms: a model with the derived average PM convolved with the random errors provided by the astrometric solution. The difference shows the amplitude of the systematic on the mean value of the PM for a bin of this size.
In the text

Fig. 15 Comparison to literature values of the proper motions of the dSph. Symbols with the same colour correspond to the same dSph, where filled symbols with error bars are those derived in this paper, and open symbols surrounded by ellipses correspond to the literature values. The error bars have the size of ${ϵ}_{{\mu }_{\alpha }*}$ and ${ϵ}_{{\mu }_{\delta }}$, as reportedin Table C.2. The black cross in the bottom left corner indicates our estimate of the systematic uncertainty on the PMs.

In the text

Fig. 16 PMs of stars in the LMC (upper) and SMC (lower) showing the components μx (left) and μy (right), which are described in the text (Eq. (3)). The colour shows the median PM in each pixel (after filtering). The black density contours are logarithmically spaced, such that the outermost contour is at a source density 100 times lower than the highest density. The centre of each colour bar is chosen to be the median PM of all sources. This is not the same as the PM derived for the clouds below, because the sources are not distributed symmetrically around the assumed dynamical centre; the photometric and dynamical (from the HI disc) centres are offset from one another. Trends in PM, particularly the trend associated with rotation in the LMC, are clearly visible. The banding associated with the Gaia scanning law, and as seen in the parallaxes, are clearly visible as well.

In the text

Fig. 17 Residual PM, after subtraction of a model PM field for the LMC (top) and the SMC (bottom). For both objects, the majority of the variation disappears, and the banding in PM is more clearly visible. Black density contours are spaced in the same way as in Fig. 16. For the LMC, the imprint of orbits on the bar can be seen as the bluer area on the lower side of the bar (as it appears in the plot) and the redder areaon the upper side of the bar. The model that is subtracted is fit from sources within angular radii ρmax = 3° of the centre for the LMC, and ρmax = 2° of the centre for the SMC.

In the text

Fig. 18 Rotation curve (top) and median vR (bottom) of the LMC. The assumed values for the centre-of-mass velocity and orientation of the disc (i and Ω) are taken from a fit to all stars within angular radii ρ < 3° of the LMC centre. Angular distances and velocities given on the lower and left axes have been converted to real-space values on the upper and right axes assuming a distance to the LMC of 50.1 kpc (Freedman et al. 2001). The points shown in the upper panel are derived from observed line-of-sight velocities of old and young stars by van der Marel & Kallivayalil (2014, their Table 4.)

In the text

Fig. 19 Examples of the orbits of some of the globular clusters in our sample. The different colours correspond to the different potentials: Model-1 (based on McMillan 2017, in blue), Model-2 (based on Allen & Santillan 1991, in black), and Model-3 (based on Robin et al. 2003, 2012, in red). The orbits of clusters that remain in the inner few kpc are quite different for the various potentials, while as expected, the differences are much smaller for those that have pericentres greater than ~ 2 kpc (e.g. NGC6496, fourth row). In these cases, the location of streams, if present, can be predicted much more reliably.

In the text

Fig. 20 Distribution of orbital parameters for the globular clusters in our sample. Three globular clusters are not shown in these plots: NGC 3201 and NGC 4590, which have apocentres ~ 30 kpc and large retrograde and prograde motions, respectively (|Lz| > 2000 kpc km s−1), and NGC 5466, which has an apocentre ~50 kpc (given its orbital parameters, it was likely associated with the Sagittarius dwarf, as proposed by Bellazzini et al. 2003). We note the concentration of clusters with small pericentre and apocentre around Lz ~ 0 for Model-3. Such clusters are on box orbits in the barred potential of Model-3, and hence do not conserve Lz, whose time average over the 10 Gyr of integration is approximately zero. The error bars correspond to the 16th and 84th quantiles obtained from the orbit integration using the statistical errors. The effect of the systematic error on the orbits of the globular clusters is found to be negligible and is not shown here.

In the text

	Fig. 21 Orbits integrated backward in time for 2.5 Gyr for the different dwarfs shown in different colours using the potential of Model-2.
In the text

Fig. 22 Distribution of orbital parameters for the dSph. The different colours indicate computations with different potentials. The agreement is generally good, with Model-3 systematically leading to smaller pericentres as a result of its higher mass at the radii probed by the systems. Larger differences are found for the more distant objects, revealing the sensitivity of their orbits to variations in the assumed mass distributions for the Galaxy. The dSph eccentricity distribution differs from that shown in Fig. 20 for the globular clusters in our sample. The symbols with solid error bars correspond to the median and uncertainties derived using the Monte Carlo realisations, and those with dotted error bars also take a systematic error of 0.035 mas yr−1 on each of the PM components into account.

In the text

Fig. 23 Top: orientation of the time-averaged orbital plane, defined by the angular momentum angles (ϕ, θ) for the different dwarfs for the various potentials. The colours and error bars are the same as in Fig. 22, and open and solid circles are the median values obtained by including or excluding our estimate of the systematic error on the PMs. Most dwarfs have highly inclined orbits with respect to the Galactic plane (i.e. θ ~ 90o). Their variation in orientation (angle ϕ) over the 10 Gyr of integration is much smaller than the size of the error bars. Bottom: histogram showing the present-day inclination of the orbits for the globular clusters.

In the text

	Fig. 24 PMsof stars in the LMC, represented as vectors, overlaid on a representation of the source density. This figure shows the clear consistent rotation measured by Gaia around the centre of the LMC.
In the text

	Fig. 25 Distribution of DR2 stars in a circle of 15 degrees radius around the Galactic anticentre direction in the L⊥ = Ly vs. Lz space. The Helmi stream is the distinct density enhancement near the bottom of the plot (in blue). We have assumed here vrot = 220 km s−1 for the LSR velocity, vz,⊙ = 5 km s−1, and vϕ,⊙ = 7 km s−1 for the peculiar motion of the Sun, and R⊙ = 8.1 kpc.
In the text

	Fig. 26 Hertzprung-Russell diagram using Gaia DR2 G magnitudes, colours, and parallaxes in the direction of the Galactic anticentre. The yellowish cloud represents a random sample of disc stars, selected on the basis of their disc-like Lz and L⊥, while the blue dots are the 32 new members of the Helmi stream.
In the text

	Fig. A.1 PM of stars within 0.5 deg. from the centre of NGC2298 as extracted from the GACS archive. Left: all ~ 259 000 stars. Middle: zoom-in on the PM field, which shows a concentration of the stars near position (3.2, −2.2) mas yr−1 corresponding to cluster members. Right: another zoom-in on the PM field for the cluster.
In the text

	Fig. A.2 From left to right: distribution of positions on the sky for the first selection on cluster star PMs, showing the cluster in the centre, PM distribution for the cluster field, and distribution of parallax vs. PM in RA.
In the text

	Fig. A.3 HR diagram for the cluster field before (left) and after (right) the astrometric solution.
In the text

Fig. A.4 Two examples of the parallax and PM dispersion diagrams for clusters at a distance of approximately 2 (top) and 5 (bottom) kpc. The parallax and PMs are plotted against the standard uncertainties. The blue lines show the expected uncertainty levels relative to the mean value at –2, –1 0 1 and 2 σ, including noise contributions from the internal velocity dispersion and the dispersion of the parallaxes from the depth of the cluster. At low levels of standard uncertainty, the contributions from the internal dispersion are easily detected. The colouring ofthe data points reflects the error correlations: red means strongly negative, blue means strongly positive, and black means low.

In the text

	Fig. A.5 Three examples of the PM dispersion in the radial and transverse directions as a function of distance to the cluster centre. The three vertical lines indicate, from left to right, the core radius, the half-light radius, and the tidal radius following Harris10.
In the text

	Fig. A.6 Two examples of astrometric data coverage with five-parameter solutions. Left: ω Cen, the worst case; right: NGC 5272, a more average example of coverage. The gaps in the coverage for ω Cen are the result of the filters that have been applied to the astrometric data. The cyan circles are at intervals of 35 pc in ω Cen and 10 pc in NGC 5272.
In the text

Fig. A.7 From left to right: standard uncertainties on the astrometric parameters of individual stars as a function of magnitude in ω Cen and NGC 5272, and standard uncertainties on the photometric data for NGC 5272. The photometric data clearly show the variability of the RR Lyrae stars (around G = 15.3), while the astrometric data also show locally much increased uncertainties at the same brightness.

In the text

	Fig. A.8 Parallax of the stars in the Carina dSph as function of their G magnitude. A clear systematic offset is apparent and is significantly measured. The expected parallax (based on literature values) is indicated by the red dashed line.
In the text

	Fig. A.9 Measured parallaxes of stars in the LMC (upper panel) and SMC (lower panel). The banding associated with the Gaia scanning law is clearly visible.
In the text

Fig. A.10 Parallaxes and residual PMs of sources in the LMC. The lines show the median value as a function of coordinate position x for all LMC sources within 0.2° stripes centred at y = −0.5, 0 or 0.5°. Residual PMs are calculated after subtraction of a model disc with parameters determined using only sources within 3° of the centre. The non-zero values are due to both the systematic errors in the Gaia data and the differences between the simple disc model and the true dynamics of the LMC.

In the text

	Fig. A.11 Correlations in the PMs of the stars in the field of view towards the different dSph galaxies in our sample. The different colours indicate the amplitude of the correlations.
In the text

	Fig. D.1 Circular velocity curves for the three Galactic potentials considered for the orbit integrations in Sect. 6.
In the text

Fig. D.2 Orbits of the globular clusters shown in Fig. 19, now plotted in cylindrical coordinates and as a function of time (right panels) for shorter integration times (approximately three radial oscillations). The different colours correspond to the axisymmetric potentials of Models 1 (blue) and 2 (black) and to the barred Model-3 (red). The solid circles in the leftmost panels denote the present-day positions, while the squares are the positions for each of the models at the end of the chosen integration time.

In the text

	Fig. D.3 Distribution of orbital parameters for the dwarfs (in different colours) showing the results obtained from the different Monte Carlo realisations (only for those within the 1σ uncertainties on the observables), computed using the potential of Model-1.
In the text