Gaia early DR3 systemic motions of Local Group dwarf galaxies and orbital properties with a massive Large Magellanic Cloud.

eGDR3 bulk motions of Local Group dwarf galaxies

G. Battaglia
S. Taibi1,2
G. F. Thomas1,2
T. K. Fritz1,2

Instituto de Astrofísica de Canarias, Calle Vía Láctea s/n, E-38206 La Laguna, Tenerife, Spain

Universidad de La Laguna, Avda. Astrofísico Fco. Sánchez, E-38205 La Laguna, Tenerife, Spain

We perform a comprehensive determination of the systemic proper motions of 74 dwarf galaxies and dwarf galaxy candidates in the Local Group based on Gaia early data release 3. The outputs of the analysis for each galaxy, including probabilities of membership, will be made publicly available. The analysis is augmented by a determination of the orbital properties of galaxies within 500 kpc.

We adopt the flexible Bayesian methodology presented by McConnachie & Venn (2020), which takes into account the location of the stars on the sky, on the colour-magnitude diagram and on the proper motion plane. We apply some modifications, in particular to the way the colour-magnitude diagram and spectroscopic information are factored in, e.g. by including stars in several evolution phases. The bulk motions are integrated in three gravitational potentials: two where the Milky Way is treated in isolation and has a mass 0.9 & 1.6 $\times 10^{12}$M$_{\odot}$ and the time-varying potential by Vasiliev et al. (2021), which includes the infall of a massive Large Magellanic Cloud (LMC).

We are able to determine bulk proper motions for 73 systems, and we consider reliable 66 of these measurements. For the first time, systemic motions are presented for galaxies out to a distance of 1.4 Mpc, in the NGC 3109 association. The inclusion of the infall of a massive LMC significantly modifies the orbital trajectories of the objects, with respect to orbit integration in static Milky Way-only potentials, and leads to 6 galaxies being likely associated to the LMC, 3 possibly associated and 1 recently captured object. We discuss the results of the orbit integration in the context of the relation of the galaxies to the system of Milky Way satellites, implications for the too-big-to-fail problem, impact on star formation histories, and tidal disruption.

§ INTRODUCTION

Knowledge of the bulk motions of galaxies residing in the Local Group (LG) is a precious
resource for a wealth of galaxy evolution and near-field cosmology investigations, e.g. inferences of the mass, barycenter position and velocity of the LG <cit.>, studies of the possible history of past interactions between the Milky Way (MW) and the M31 system and its future fate <cit.>; determinations of the mass of the MW through dynamical modelling of tracers of its gravitational potential, such as its satellite galaxies <cit.>; group infall as well as the significance and stability of the Vast Polar Structure <cit.>; considerations on the missing satellite problem <cit.>. The orbital history of MW satellite galaxies is also very likely to influence several aspects of their evolution, through e.g. the impact of ram-pressure stripping and tidal effects onto their gas content, star formation history (SFH), morphology and dark matter (DM) halo properties <cit.>.

Before the second data release of the Gaia mission (GDR2) [Gaia Collaboration et al., 2016, Gaia Collaboration et al., 2018], measurements of the systemic proper motions (PMs) of galaxies in the LG were essentially limited to the Magellanic Clouds, the so-called classical MW dwarf spheroidal galaxies (dSphs), one "ultra faint dwarf" (UFD), M31, M33 and IC 10, mostly from HST observations and a few VLBI observations of OH masers (see references in Sect. <ref>); some of them having become available surprisingly recently <cit.>.

Since April 2018, the situation has seen a dramatic improvement, starting with the Gaia science verification article [Gaia Collaboration et al., 2018]. Multiple determinations of the systemic PM of a large number of MW satellite galaxies and galaxy candidates blossomed in a matter of weeks after GDR2, led by several groups in the community, and using a variety of techniques, e.g. focusing only on stars with prior spectroscopic information [Simon, 2018, Fritz et al., 2018] or including the full set of stars with astrometric information [Kallivayalil et al., 2018, Massari & Helmi, 2018]. It is now becoming routine to use Gaia astrometric data also to remove contaminants, as well as to attempt systemic proper motions determinations along with the study of other properties of the systems <cit.>. Surveys of the MW stellar halo and sub-structures within are and will make plentiful use of Gaia astrometry to boost the success rate in target selection <cit.>.

The methodologies applied in the early GDR2 works mentioned above were rather simple ones, being based on iterative cleanings of the data-sets via $\sigma$-clipping and no statistical treatment of the foreground/background contamination. Later on, more sophisticated methods were used, e.g. with simultaneous statistical modelling of the properties of the dwarf galaxy and the contamination. For instance, Pace & Li, 2019 used the spatial and PM  information of stars preselected to have magnitude and color lying on an isochrone and adopted a multi-variate Gaussian in proper motion for both the dwarf galaxy and the MW, while McConnachie & Venn, 2020 used all the  observables at once and adopted the empirical distribution of the contaminant stars in the PM and the colour-magnitude (CM) planes.

Interestingly, the use of GDR2 data has been also pushed beyond the MW system, with determinations of the tangential motions of M31 & M33 [van der Marel et al., 2019], as well as of a few LG dwarf galaxies such as NGC 6822, IC 1613, WLM and Leo A  [McConnachie et al., 2021].

The early third data release of Gaia data, hereafter eGDR3,  [Gaia Collaboration et al., 2021] has implied more precise and accurate astrometric measurements; in particular, for PMs the precision has increased of a  factor of two and systematic errors decreased by a factor $\sim$2.5 [Lindegren et al., 2020]. McConnachie & Venn, 2020 provides updated systemic PMs for the 58 MW satellites previously considered by the same team with DR2 data, with the improved astrometry now allowing to detect the systemic PM of Boötes IV, Cetus III, Pegasus III and Virgo I. Recently, Li et al., 2021 provide an independent determination of systemic PMs for 46 MW satellites, and integrate the 3D motions in four isolated MW potential models, with a  total mass from $2.8 \times 10^{11}$ to $15 \times 10^{11}$ M$_{\odot}$. Martínez-García et al., 2021 combined the astrometric and spectroscopic information available for 14 MW satellites to study their internal kinematics and quantify the presence of velocity gradients.

In this work we aim at providing a comprehensive determination of systemic proper motions based on eGDR3 not only for MW satellites, but for LG dwarf galaxies in general[For brevity, hereafter we refer to all the systems in the sample as "dwarf galaxies", including when they are only "candidates" and even in the case of larger galaxies like M 33 and NGC 3109.], and to push for the first time these determinations to even larger distances, i.e. reaching out to the NGC 3109 association, at $\sim$ 1.4 Mpc. We are making use of the best techniques in the literature, i.e. those by McConnachie & Venn, 2020, McConnachie et al., 2021, inspired by Pace & Li, 2019, since they take full advantage of the observables available for the largest number of stars with full astrometric and photometric information (location on the sky with respect to the dwarf galaxy centre, and on the colour-magnitude diagram and on the proper motion planes), and model  them in a Bayesian way with a mixture model accounting for contaminant sources. We have introduced a few modifications to these techniques, mainly aimed at making an even more realistic treatment of the stellar population content of the dwarf galaxies and its distribution on the colour-magnitude plane accounting also for the photometric completeness of eGDR3 data; this allows us also to determine probability of memberships for stars in different evolutionary phases, which we will make available to the community [http://research.iac.es/proyecto/GaiaDR3LocalGroup/pages/en/home.php], together with several other outputs of our analysis.

Finally, we study the orbital properties of the galaxies surroundings of the MW by integrating their bulk motions in a set of gravitational potentials, bracketing a 0.9-1.6$\times 10^{12}$ M$_{\odot}$ range for the mass of the MW. Motivated by the work by Patel et al., 2020, who showed that the orbits of MW satellites can differ significantly between a gravitational potential including only the MW and one where the gravitational influence of the LMC (and SMC) are taken into account (and all galaxies are free to move in response), we integrate the bulk motions also in the triaxial time-varying MW potential made available by Vasiliev et al., 2021, where the infall of a massive LMC and the response of the MW to this infall are modelled. In this context, we also revisit the association of dwarfs surrounding the MW to the LMC system.

In Sect. <ref> we introduce the sample of galaxies analysed and in Sect. <ref> the data-sets used and the quality selection criteria applied. In Sect. <ref> we present the methodology for the systemic proper motion determinations, and discuss the output and the tests performed to tackle the robustness of the method; in Sect. <ref> we complement the results with a determination of the zero-points and additional errors due to systematics in the Gaia eDR3 data, using QSO. Our systemic proper motions are compared to those in the literature in Sect. <ref>. In Sect. <ref> the bulk motions are integrated in the three gravitational potentials and the resulting orbital trajectories and parameters are then used to address the impact of the LMC on the reconstructed orbital history, make considerations on the too-big-to-fail problem, on the system of MW and LMC satellites, and on observed properties such as SFHs.  We discuss other potential applications of our work in Sect. <ref> and present our conclusions and summary in Sect. <ref>.

§ SAMPLE

L.o.s. velocity as a function of heliocentric distance for the galaxies in the sample (circles). The systems without a literature measurement of the l.o.s. velocity are assigned a null value, exclusively for the purpose of this figure, and are indicated as crosses. The color-coding is based on the uncertainty in heliocentric transverse velocity, derived from the statistical uncertainties in the systemic proper motions and distance modulus (see Sect. <ref>).

The 74 systems considered in this work are listed in Tab. <ref>, together with their main global properties (see also Fig. <ref>). The sample is the union of the dwarf galaxies studied in Fritz et al., 2018, Fritz et al., 2019, and some other recently discovered satellites of the MW. We also included isolated LG dwarf galaxies within $\sim$1.4$\,$Mpc.

Due to the large distance of the isolated LG dwarf galaxies, in terms of resolved or partially resolved sources, only H II regions, young main sequence, blue/red super-giants and asymptotic giant branch (AGB) stars brighter than the tip of the red giant branch (RGB) are potentially detected above the magnitude limit for Gaia astrometric measurements. Therefore not all isolated dwarf galaxies will have enough (or any) eGDR3 sources  with full astrometric solution.

In order to select which systems to consider on a first pass, we gave priority to those galaxies within 1.5Mpc with an HI detection, as listed in McConnachie, 2012,  and hence likely to host young stars and/or HII regions. The systems that showed to have a clear enough detection in Gaia eDR3 of centrally concentrated sources (after a first, rough, spatial, parallax and proper motion selection) were retained. We exclude IC 10 because, apart from being only marginally detected, its analysis is complicated by the very high and patchy extinction  in its direction. The borderline systems like LGS3, Antlia, NGC 205, NGC 185 were excluded but the closer galaxy Leo T was retained.

We exclude the Magellanic Clouds, M 31 and the Sagittarius (Sgr) dwarf galaxy, because we are neglecting internal motions in our analysis, while they might be relevant for the systemic PM  determination of these systems. These galaxies have anyway already been the subject of very detailed analyses based on Gaia data <cit.>. Note that in their GDR2-based study, van der Marel et al., 2019 found no difference in the centre-of-mass proper motion when explicitly modelling M33 rotation or neglecting it. This means we include the following galaxies outside the virial radius of the MW: Eridanus II, Leo T, Phoenix, NGC 6822, WLM, IC 1613, Leo A, M33, Peg-dIrr, UGC 4879, Sgr-dIrr, Sextans A, Sextans B, NGC 3109.

We caution the reader that, even for those galaxies for which a determination of the systemic PM is possible, depending on the distance, the error in the transverse velocity is still too large for scientific applications (see Fig. <ref>). We refer the reader to Sect. <ref> to learn about which galaxies have their error budget dominated by random or statistical errors.

§ DATA-SETS

Gaia eDR3 astrometric and photometric measurements for stars with full astrometric solutions constitute the bulk of the data over which the analysis is performed.
For the distant systems (i.e., those beyond 400 kpc) we complement Gaia photometry with the deeper data-set from the Pan-STARRS1 Surveys <cit.>; this will be used to identify the regions in the colour-magnitude plane where to select candidate massive blue stars and red super giants stars and to determine their spatial distribution, in some cases.
We refer to Sect. <ref> for a detailed explanation of the methodology.

We download eGDR3 and PS1 data over an area centered on the systems under consideration, with a radius such to guarantee at least 2000 objects where we expect contaminants to be clearly dominant over stars from the galaxy under consideration, i.e. beyond 5 half-light radii, $R_h$[The colour-magnitude diagrams were inspected for all the systems, to confirm the dominance by contaminant sources in those regions.].

We concentrate exclusively on those eGDR3 detected objects that are not flagged as duplicated source, have a full astrometric solution (astrometric_params_solved $\ge$ 31),  high-quality  astrometry  <cit.> and reliable photometric measurements. For this last factor, we retain the measurements with an absolute value of the corrected excess factor within 5 times the standard deviation expected at the corresponding G-mag   <cit.>. In order to exclude sources seen as extended or not isolated by Gaia, we retain objects with ipd_frac_multi_peak $\le$ 2 and   ipd_gof_harmonic_amplitude $<$ 0.2; this latter cut is less restrictive with respect to what adopted in Fabricius et al., 2021 but adjusts better to the distribution of values seen for the sources around the  galaxies under consideration. Objects with source_id with a match in the Gaia AGN catalogue are excluded.

We exclude clear foreground stars by requiring the parallax of each individual source to be consistent with the parallax expected at the dwarf distance modulus within 3$\sigma_\pi$; here $\sigma_\pi$ is the sum in quadrature of the parallax error on the individual measurements and that due to the uncertainty of the galaxy distance modulus. We apply a global zero-point offset of $-0.017$ mas to the parallax measurements of the individual stars [Lindegren et al., 2020]. We do not correct for the
Gaia parallax zero-point as a function of location, magnitude and colour, because of its negligible effect on our analysis: even at the brightest magnitudes considered here ($m_G \sim 13.4$ mag for the tip of the RGB of Delve 1, the closest system in the sample), the maximum difference between the zero-point applied and that expected at $m_G \sim 13.4$ mag would be $\sim$0.03 mas [Lindegren et al., 2020], smaller than the $3\times \sigma_\pi$ range under consideration.

The apparent $G$-mag for the sources with 6-parameters solutions are corrected as in Riello et al., 2021, with the python code presented in Gaia Collaboration et al., 2021. Finally, the apparent $G-$mag and the $BP-RP$ color are corrected for extinction using the Schlegel et al., 2011 maps interpolated at the position of the stars and using the  Marigo et al., 2008 coefficients for the Gaia filters, based on <cit.> <cit.>.

§ DETERMINATION OF SYSTEMIC PROPER MOTIONS

§.§ Method

Rather than relying only on eGDR3 sources with previous spectroscopic observations, we adopt the flexible methology by <cit.>; <cit.>, inspired by Pace & Li, 2019,
which allows to make use of all the stars with Gaia astrometric and photometric measurements (aside from the quality cuts detailed in the previous section). We refer the reader to the original source for a detailed explanation of the methodology. Here suffices to say that it is based on a maximum likelihood procedure, which has three  free-parameters: the systemic proper motion of the galaxy $\mu_{\alpha,*, sys}$ and $\mu_{\delta, sys}$ and the fraction of stars in the dwarf galaxy under consideration over the total, $f_{gal}$. It is assumed that  the intrinsic dispersion is negligible in the distribution of PM measurements.

The likelihood of a star to belong to the system/object  under consideration, $L_{\rm gal}$,  or to the contamination, $L_{\rm c}$, is estimated taking into account a spatial, colour-magnitude and proper motion likelihood term.

After having determined $\mu_{\alpha,*, sys}$, $\mu_{\delta, sys}$ and $f_{gal}$, these, together with the likelihoods, can be used to obtain the probability of membership of each star to a given galaxy (Eq. 5 in MV20a):

\begin{equation}
P_{\rm gal} = \frac{f_{\rm gal} L_{\rm gal}}{f_{\rm gal} L_{\rm gal} + (1 - f_{\rm gal}) L_{\rm c}}
\end{equation}

Below we describe the methodology followed to determine the various terms of the likelihood function, and provide an overview of how each galaxy was treated in Tab. <ref>.

§.§.§ Spatial Distribution

The contaminants are assumed to be uniformly distributed over the areas considered around each system.

The spatial term of the likelihood function for the stars belonging to the dwarf galaxy is based on the 2D structure of the dwarf galaxy stellar component. In order to evaluate it, we adopt two main approaches, depending on the system under consideration: either we parametrize it to have an elliptical shape and an exponentially declining surface number density profile, as done by MV20a, or determine it empirically, as done by Mc21 ("Exp" and "Emp" in column "Spatial" of Tab. <ref>).

We refer to MV20a and Mc21 for the possible caveats concerning this approach.

For the "Exp" case, a 2D look-up map is created by co-adding and then normalising $N=1000$ MonteCarlo realizations of the expected 2D surface number density at a given position on the sky; in each realization, values for the ellipticity, position angle and half-light radii[The ellipticity is defined as $1 - b/a$, where $b$ and $a$ are the projected minor and major axes; the position angle increases from North to East; the half-light radius here refers to the projected one on the sky, along the major axis.] are randomly extracted from a Gaussian distribution centred on the values listed in Tab. <ref> and with dispersion given by the average of the lower and upper 1-$\sigma$ uncertainties listed in the same table. For those galaxies where a determination of the ellipticity is missing or only an upper limit is available, we assume the spatial distribution to be circular. This approach is applied essentially to all the systems for which the spatial distribution of the majority of sources detected by Gaia and that we will use in the analysis have a smooth, spheroidal-looking morphology (i.e. the "classical" dSphs and the UFDs) or for some late-type systems for which there is not enough statistics for an empirical determination (see below).

For the distant late-type systems, the majority of the sources detected by Gaia will be blue massive stars, red super-giants (RSGs) and AGB stars. Among these, we will concentrate on the blue and RSGs as more easily identified on the CMD (see also Mc21). These young stars are those that give an "irregular" morphology to some of these galaxies, due to asymmetries in their spatial distribution. Therefore, the approach to follow for the late-type systems is decided after a visual analysis of the spatial distribution of stars with colors consistent with being young main-sequence or blue super-giant stars in the PS1 photometry: if their spatial distribution is well defined, then the probability distribution of the spatial term is determined empirically from these stars as a normalized 2D histogram within 3 half-light radii (apart from M33, that is missing this quantity, for which we consider 0.5 deg). For those cases where the statistics of blue stars are not sufficient for an empirical determination, we resort to modeling the spatial term as an exponentially declining profile. In all cases, the spatial distribution of the RSGs is assumed to follow that of the blue stars, which is true to a good approximation.

§.§.§ Distribution on the colour-magnitude plane

Spatial     CM     Galaxies

Exp      Emp     Carina, Draco, Fornax, Leo I, Leo II,

Sculptor, Sextans, Ursa Minor

Exp     Syn      All those not listed in the other rows

Emp     Box     IC 1613, NGC 6822, WLM

NGC 3109, Sextans A, Sextans B, M33

Exp     Box     Leo A, Peg-dIrr, Sg-dIrr, UGC 4879

Methodology used for the spatial and CM  term of the likelihood for stars in the dwarf galaxy, when determining the systemic PM (Sect. <ref>. In the column "Spatial", "Exp" indicates that an exponential profile was used for the surface number density distribution and "Emp" that the spatial distribution was determined empirically. In the column, "CM", "Emp" indicates that the probability distribution on the CM plane was determined empirically (from within 1 half-light radius), "Syn" that the synthetic CMD was used and "Box" that a uniform probability was given to stars within regions of the CMD compatible with hosting massive blue stars and red super-giants.

We concentrate on sources with $-1.0<$(BP-RP)$_0 <2.5$,
apart for the late-types, for which we adopt $-1.5<$(BP-RP)$_0 < 2.5$.
We have verified that these color cuts works also for metal-rich systems such as Fornax, WLM, NGC 6822 etc.

As in MV20a, the distribution of contaminants onto the colour-magnitude plane is determined empirically, from eGDR3 sources at semi-major axis radii larger than $5\times R_{h}$, with the exception of those systems for which this limit exceeded the spatial extension of the catalogue downloaded (e.g. for Antlia II), in which case we assume radii larger than $3\times R_{h}$.
For M33 we do not use a value of the half-light radius but define the region by eye, beyond 1deg.

For the CM probability distribution of stars belonging to the systems of interest, we introduce a few changes with respect to the method of MV20a:

* Empirical determination ("Emp" in column "CM" of Tab. <ref>) For well populated and nearby systems, such as the "classical dSphs", the distribution on the CM-plane is determined empirically, from the region within one half-light radii. The dwarf galaxy's stellar population within this region dominates over that of the contaminants and the increase in statistics over considering a smaller area is worth the introduction of a few contaminants.  While we do not expect the choice of using an empirical determination of the dwarf's CM likelihood term to cause a significant difference in the determination of the systemic proper motion over e.g. using an isochrone, we wish to factor in the CMD information in the estimate of the probability of membership for stars in different evolutionary phases, since this is one of the products that we make available.

Classical dSphs are well-known to display stellar population gradients; however a complete modelling of the CMD as a function of distance from the dwarf centre is outside of the scope of this work, and the CMD of the central regions contains all the features present also in the outer parts.

* Synthetic CMD ("Syn" in column "CM" of Tab. <ref>) All those galaxies closer than 440 kpc that are not included in the category above do not have enough signal inside their half-light radius for an empirical determination of the CM probability distribution. In addition, many of these systems are very faint and sparsely populated;  therefore we wish to adopt an approach that includes all the relevant evolutionary phases and does not exclude a priori possible members, e.g. if they were not to fall on the locus of an isochrone of a given age and metallicity.

To this aim, rather than an isochrone as done in MV20a, we use Basti-IAC to create a synthetic CMD in the eGDR3 photometric filters[The synthetic CMD was courtesy of S. Cassisi; the website of Basti-IAC is  http://basti-iac.oa-abruzzo.inaf.it/index.html.], based on the stellar evolutionary models presented in Hidalgo et al., 2018. These models include also the He-burning phase, which can be precious in faint systems, since in particular the blue part of the horizontal branch is a region of rather reduced contamination.  An advantage of using a synthetic CMD over using an isochrone is that stars are distributed on the magnitude and color plane in the correct proportion  (for a given SFH, chemical enrichment law, initial mass function,..).

Given that the systems we are applying this method to are either completely or mostly dominated by ancient stellar populations, we adopt a constant star formation rate between 12 and 13 Gyr ago, and a metallicity centered around [Fe/H] $=-2.3$ with a spread of 0.5 dex; this is representative of the metallicity distribution function of stars in UFDs <cit.>.

One hundred realizations are carried out, where the synthetic CMD is shifted in distance modulus, drawing from Gaussian distributions centered on the values listed in Tab. <ref> and with dispersion given by the average of the upper and lower 1-$\sigma$  errors. At the same time, the photometric errors are introduced by scattering the BP-RP colours of the stars in the synthetic CMD according to the photometric errors derived from the eGDR3 catalogue corresponding to each given object, at the appropriate apparent G-mag.

Correction for photometric completeness The synthetic CMD of course does not suffer from photometric incompleteness; on the other hand, it should be considered that the CMD of the contaminants does suffer from this issue, since it is derived empirically from the eGDR3 data, and that the completeness varies depending on the number of transits in a given region of the sky. This might alter the relative probabilities of dwarf galaxy's stars versus contaminant stars in some regions of the CMD, in favour of the former.

In order to introduce an (approximate) corretion to take this effect into account, we resort to Gaia Universe Model Snapshot <cit.>. We download GUMS models from the Gaia archive around the position of each system, and calculate the ratio of the luminosity function of contaminants stars in the observed eGDR3 catalogues (which will be mostly MW stars) and that of model stars, in a representative color range, 0.5 $<$ (BP-RP)$_0$ $<$ 2.5. In the assumption that the model is a reasonable approximation of the data, at G-mag where Gaia should not suffer from completeness issues, this ratio should be around unity. In practise, there are some deviations; therefore, we normalize the ratio to the median value in the range 17 $<$ G-mag $<$ 20. After this step, the ratio of the two luminosity functions oscillates around one, except at faint magnitudes, where a decline towards zero is observed, which can be assumed to be due to the Gaia photometric completeness. This will be the factor by which we multiply the counts in the CMD look-up map as a function of magnitude[In practise, to avoid introducing noise, we only multiply the counts by this correction factor at G-mag $>$19 and if the ratio is $<$0.75.]. As we will see in Sect. <ref>, this correction has a minor effect, but this will be adopted for our baseline results, listed in  Tab <ref>.

Even though Phoenix and Leo T do contain a sprinkle of young, blue stars detected in eGDR3, since the majoritarian population in eGDR3 data is by far represented by RGB stars, these systems are treated with the synthetic CMD. This choice does not impact the determination of the systemic PM but it implies that these young stars will be missing from our list of probable members.

* Box ("Box" in column "CM" of Tab. <ref>) For the distant ($>$ 440 kpc) and well populated galaxies, we follow closely the approach by Mc21, and focus on blue sources and RSGs. The color and magnitude limits are chosen by visual inspection of the PS1 photometry, and transferred to the Gaia eGDR3 bands, using the stars in common between the data-sets for each galaxy. We assign a uniform probability inside these boxes.

For both the dwarfs and the contaminants, we construct the CM  look-up map in bins of magnitude and color, smooth it with a boxcar kernel and then proceed to normalizing it.

§.§.§ Proper motion

For the distribution on the PM plane, we adopt the same approach as MV20a, i.e. an empirical determination for the contamination, while assuming a multi-variate Gaussian distribution for the dwarf galaxy, taking into account the correlation terms between the $\mu_\alpha,*$ and $\mu_\delta$ of the individual stars.

For the MW classical dSphs and UFDs, we restrict the range of the analysis to within $\pm$5 mas yr$^{-1}$, corresponding to a generous tangential velocity cut of $>$470 km s$^{-1}$ at a heliocentric distance larger than 20 kpc ($\pm$3 mas yr$^{-1}$ for Antlia II, to reduce the overwhelmingly large contamination).

For the other systems we filter out sources whose proper motion in each component at the distance of the galaxy would imply tangential velocities 3x in excess of a given velocity dispersion (around the reflex proper motion at the galaxy's sky location). As dispersion, we consider the square-root of the quadratic sum of the uncertainty given by the proper motion measurements and 200 km s$^{-1}$, where the latter is the observed scatter in l.o.s. velocities for the whole sample of galaxies in Tab. <ref> (assuming that the scatter in tangential velocities is the same)[We tested the performance when relaxing the cut in tangential velocity and the resulting systemic PMs are always within 1 or at most 2-$\sigma$ from each other, showing that the method is robust also when removing this condition.].

§.§ Results

The outcome of the analysis is summarized in Tab. <ref>.

For the majority of the systems in the sample, in output there are sizeable numbers of stars with large membership probability (42 and 52 systems with $>$10 stars with P$>0.95$ and P$>0.5$, respectively). In fact, some galaxies are extremely well populated. For example, we obtain $>$1000 stars with membership probabilities P$>$0.95 in each of the classical MW dSphs (except Leo II) as well as in M 33 ($>$20000 for Fornax). Ten systems have between 100-1000 $P>0.95$ members, among which several of the distant galaxies (Phoenix, NGC 6822, IC 1613, WLM, NGC 3109).

On the other end of the spectrum, there are systems with only an handful of probable members or none at all. The analysis does not lead to a systemic PM determination for Pisces II and Virgo I, with PDFs that are essentially flat. Other systems with clearly problematic PDFs are DESJ0225+0304, Pegasus III, Tucana V, with strong lopsideness and/or very extended wings of high amplitude, and Indus I, with a double peaked PDF. For Pisces II and Tucana V, however, we are able to obtain a systemic PM when including the information about l.o.s. velocities for the stars observed spectroscopically (see Sect. <ref>). We would advice against using the systemic PMs for all these cases (whose names are highlighted in red in Tab. <ref>) and advice to use the motions obtained when considering the spectroscopic information for Pisces II and Tucana V.

In addition to the above, the shape of the PDFs and of the distribution of probable member stars on the plane of the sky, proper motion and colour-magnitude leads us to advise exerting caution when considering the motions of Cetus III, Indus II, Aquarius II, Delve 1, Reticulum III, Bootes IV (see Sect. <ref> for more details; these are highlighted in orange in Tab.<ref>).

Example of plots of the distribution of member stars projected on the tangent plane passing through the galaxy centre, and on the proper motion and colour-magnitude plane are given in Figs.<ref>-<ref>, for systems in various regimes, in terms of number of member stars, heliocentric distances, and morphological types.

In summary, out of 74 systems analysed, we are able to determine systemic proper motions for 72 systems without considering complementary spectroscopic information (73 when considering the spectroscopic information), and consider certainly reliable 62 (64 with spectroscopy). The majority of these 64 systems are found in the vicinity of the MW, within 300 kpc, but the galaxies for which we provide systemic PMs are as distant as 1.4 Mpc (Sextans A and Sextans B). This is the largest, and most extended in volume, set of systemic proper motions for galaxies and galaxy candidates. Of course, it should be kept in mind that the same uncertainty in proper motion will translate into an uncertainty in transverse velocity 10x larger for a galaxy at 1 Mpc than for one at 100 kpc!

Distribution of member stars (large circles) projected on the tangent plane passing through the galaxy centre (left), and on the proper motion (middle) and colour-magnitude plane (right), for systems in the regime of $>$500 stars with P$>0.95$, in increasing order of distance from top to bottom. The galaxy names are indicated in the figure titles. The colour-coding indicates the probability of membership (only when above $>0.5$; the stars with P $<0.5$ are shown as grey dots. The ellipses in the left panel have semi-major axes equal to 1x and 3x the half-light radii in Tab. <ref> (apart from M33, that is missing this quantity, for which we consider 0.5 deg), and ellipticity and position angle taken from the same table.

As in Fig. <ref> but for the regime of 50-200 stars with P$>0.95$.

As in Fig. <ref> but for the regime of $<$50 stars with P$>0.95$. Cetus III and Bootes IV are two cases in which the uncertainties in McConnachie & Venn, 2020 are much smaller than in our determination.

§.§.§ Tests and validations

$\bullet$ Photometric completeness As discussed in Sect. <ref>, a correction for the photometric completeness of eGDR3 data was applied to the CM probability distribution of the galaxies treated with the synthetic CMD (see Tab. <ref>). Fig. <ref> shows the comparison of the systemic proper motions determined with and without applying this correction: the determinations are always in very good agreement, well within the 1-$\sigma$ errors, apart from Antlia II, for which a larger difference is seen, but still within 2-$\sigma$. The good agreement between the two determinations is likely due to the correction kicking in at faint magnitudes, there were the individual proper motions are less accurate and therefore have a lower weight in the global determination.

$\bullet$ Inclusion of spectroscopy information MV20a analyzed the galaxies in their sample with and without considering additional information on the stars's membership from spectroscopic observations. Specifically, when including information from spectroscopy, they modified the prior on the systemic proper motion by multiplying it by a bi-variate Gaussian with mean and dispersion given by the weighted mean proper motion and associated uncertainty of the stars with Gaia astrometric measurements that are also probable spectroscopic members. The authors concluded that information from stars with spectroscopic follow-up was not required to obtain reasonable estimates of the systemic PMs; in fact, in only a few systems, the inclusion of this information played an important role: Carina III, Segue 1,  Triangulum II  and  Tucana IV, for which the spectroscopy allowed to go from a bi-modal to a uni-modal PDF (which in the case of Carina III was due to the presence of Carina II in the field-of-view). In general, the inclusion of this prior had the effect of reducing the size of the error-bars, even though in some cases only slightly. The authors also warn about the dangers of including this information for systems with only an handful of spectroscopic members, since interlopers could of course lurk  among them too.

We explore the possible improvements due to the use of spectroscopic information in a different way: rather than using the iteratively derived mean proper motion of spectroscopic member stars to modify the prior, we introduce another term in the likelihood, in which the l.o.s. velocity distribution is modelled as the sum of two 1D Gaussians, one for the dwarf galaxy and one for the MW. Since not all the stars with astrometric information do have a l.o.s. velocity measurement, we assign a l.o.s. velocity equal to 0 km s$^{-1}$ and a l.o.s. velocity uncertainty of 10000 km s$^{-1}$ to the stars in the astrometric sample that do not have a spectroscopic match; these arbitrarily large velocity uncertainties have the effect of giving these stars no weight in the estimate, but do allow us to treat the spectroscopic information as a further likelihood term.  The value of the peak l.o.s. velocity and l.o.s. velocity dispersion for the dwarf galaxies are fixed to the values in Tab. <ref>, fixing the velocity dispersion to 5 km s$^{-1}$ when only an upper limit is available.
At the same time we solve for the peak l.o.s. velocity and l.o.s. velocity dispersion of the MW component (we have tried also keeping them fixed to 0 km s$^{-1}$ and 200 km s$^{-1}$, respectively, and the results do not vary).
Fig. <ref> shows the comparison of our baseline case with the determinations using a spectroscopic prior for the category for which the largest differences could in principle be expected, i.e. the ultra faint dwarfs.
The only two systems with significant differences are Pisces II and Tucana V, where the size of the error-bars reduces  drastically. This is due to two stars with high probability of membership found for Pisces II and a PDF with wings of much lower level for Tucana V in the run with spectroscopic information with respect to that without. For the great majority of the other cases, the differences are minor both in terms of systemic motions and associated uncertainties. In Appendix <ref>, we comment on those cases where the difference between the systemic PMs with and without the spectroscopic information is larger than 0.5 $\sigma$.

As MV20a concluded, it is re-assuring that spectroscopic follow-up is not a necessary condition for systemic PM determinations.
Overall we see a lower degree of improvement than that found by  MV20a on GDR2 data. Likely, the main reason is that eGDR3 data, in particular the PMs, have become more precise, which makes it easier for the algorithm to find galaxies even if they have only a few stars above the Gaia magnitude limit. This might make it potentially easier to find galaxies using only Gaia data in the future data releases, but perhaps also already in eGDR3 <cit.>.

$\bullet$ RRLyrae As an additional check of the robustness of the results, we compared the systemic PMs with the individual measurements for RRLyrae found at projected distances within 5$\times$ half-light radii and with magnitudes approximately compatible with the horizontal branch of each system. As catalogue of RRLyrae stars, we use the union of the Gaia DR2 SOS gaiadr2.vari_rrlyrae [Holl et al., 2018, Gaia Collaboration et al., 2019, Clementini et al., 2019], the stars classified as RRLyrae of ab, c, d type in the general variability catalogues gaiadr2.vari_classifier_result, and the PS1 RRLyrae of ab or cd type by [Sesar et al., 2017][With classification score above 0.6.]. The comparison is very good. The only system in which some outliers in the RRLyrae PMs are found is Hydrus I, but this is not a cause for concern, because a closest examination shows the presence of RRLyrae compatible with belonging to the SMC-LMC system in the background and we made no attempt of statistically account for contamination in the RRLyrae variables data-set.

§ SYSTEMATIC ERRORS AND DISTANCE ERRORS

It is known that Gaia PM measurements are affected by systematic errors, which can be thought of as a component on small angular scales, $\lesssim$ 1 deg, and a component on large scales, with a scale-length of $\sim$16 deg for eGDR3 <cit.>.  Given the spatial scales typically involved in our analysis, the effect of the large-scale component should be to act as a zero-point in the observed systemic motions. On the other hand, the small-scale component will average out for some of the systems with the largest angular size, but not for a significant number of them.

Therefore we follow two routes: we treat the bias on small-scales as an additional source of noise, while we determine the zero-point from the large-scale component for each galaxy separately from QSOs.

For the small-scale error, we use the determination by Vasiliev & Baumgardt, 2021[Their equation 2, but with 400/(1+$\theta$/3), not 400 x (1+$\theta$/3), which was a typo (Vasiliev, private communication).], rather than that by Lindegren et al., 2020, since the former was derived on globular cluster stars, in which there are more close neighbours on smaller scales than among quasars, used by the latter work. As typical scale of our systems, we use the "circularized"  half-light radius, $\theta_{\rm half}$ (for M33 we use the radius containing half of the member stars). This leads to errors between 13 and 23 $\mu$as yr$^{_1}$ for both dimensions, $\sigma_{\rm vas} (\theta_{\rm half})$.

For the other component, we calculate the weighted average of the PMs of QSOs (from the table agn_cross_id provided within Gaia eDR3) within 7 deg around each galaxy. We found this scale to be a good choice in terms of overall error and scatter among the galaxies, among the explored scales of 3-10deg with 1deg steps. We concentrate on QSOs with 5p solutions, since they are known to have more precise measurements [Lindegren et al., 2020, Fabricius et al., 2021], and retain those with $G < 19$, to reduce statistical errors, and ruwe $<$1.4,  ipd_gof_harmonic_amplitude $<$ 0.2 and  ipd_frac_multi_peak $\le$ 2 for ensuring good astrometric measurements.  The zero-point, to be subtracted to the systemic PM, is calculated as a weighted mean (and its error as error in the weighted mean) after two iterations. This yields a minimum of 50 QSO, with the median being $\sim$900.

Since the Vasiliev & Baumgardt, 2021 formula includes the effect from the large-scale component, we subtract from $\sigma_{\rm vas} (\theta_{\rm half}$) the corresponding value from the same formula on the scales of the determination from QSO and account for the error on the weighted mean of QSO PMs.

Both the zero-points and additional error per PM component are given in Tab. <ref> and are used for the orbit integration analysis in Sect. <ref>) (the PMs and errors in Tab. <ref> do not include these additional errors/corrections).

In general, we find that the dominant[Here defined as being at least 1.2$\times$ larger.] source of error is the random one for all the galaxies, apart from
Fornax, Sculptor, Ursa Minor, Draco, Carina, NGC 6822, Leo I, Sextans, Antlia II, Bootes I, Hydrus I, Reticulum II, Carina II, IC 1613, Crater II, M33. Since Lindegren et al., 2020 find that the systematic PM error decreases with a similar factor with time as the random error, this is not expected to change.

In the great majority of applications, systemic PMs need to be converted into a velocity, and uncertainties in the distance modulus will contribute to the uncertainties in the physical transverse velocity. Therefore, it is interesting to known in which cases that is the largest source of error (see also Tab. <ref>); these are: Bootes I, Bootes II, Bootes III, Carina, Carina II, Carina III, Cetus II, Coma Berenices, Delve 1, Fornax, Grus II, Horologium I, Hydrus I, Phoenix II, Pictor II, Reticulum II, Sagittarius II, Segue 1, Segue 2, Tucana II, Tucana III, Tucana IV, Tucana V, Ursa Major I, Ursa Major II, Willman 1. It turns out that the majority of galaxies within 100kpc have their uncertainty in transverse velocity dominated by distance errors, as compared to that due to the statistical and systematic errors in the systemic PM.
Note that we are including an additional 0.1mag error in the distance modulus of galaxies whose published uncertainties are lower than that value; this in order to mimic the typical mismatch between values of distance modulus found from different techniques. If we drop this additional factor, the situation changes only for Fornax and Bootes II, which become dominated by the PM component. Sticking to the published uncertainties, there are systems where the distance factor can be as large as 3 to 7 times the systemic PM one (going from 10-40 km s$^{-1}$ the former, while the latter is within 3-12 km s$^{-1}$), like Carina III, Cetus II, Hydrus I, Reticulum II, Sagittarius II, Segue 1, Tucana II, Ursa Major II. These are all very faint systems, and it will be hard to improve on their distance estimates, but it might be worth the trouble.

§ COMPARISON WITH THE LITERATURE

In this section, we compare our systemic PMs determinations with those in the literature. These were obtained with Gaia DR2 data [Gaia Collaboration et al., 2018, Simon, 2018, Simon et al., 2020, Fritz et al., 2018, Fritz et al., 2019, Carlin & Sand, 2018, Massari & Helmi, 2018, Kallivayalil et al., 2018, Pace & Li, 2019, Pace et al., 2020, Carlin & Sand, 2018, McConnachie & Venn, 2020, Longeard et al., 2018, Longeard et al., 2020, Torrealba et al., 2019, Mau et al., 2020, Cerny et al., 2020, Chakrabarti et al., 2019, Gregory et al., 2020, Mutlu-Pakdil et al., 2019], eGDR3 [McConnachie & Venn, 2020, Vasiliev & Baumgardt, 2021, Jenkins et al., 2020, Martínez-García et al., 2021, Li et al., 2021], with HST [Piatek et al., 2003, Piatek et al., 2005, Piatek et al., 2006, Piatek et al., 2007, Pryor et al., 2015, Piatek et al., 2016, Sohn et al., 2013, Sohn et al., 2017] and VLBI [Brunthaler et al., 2005]. For the HST measurements, we ignore older determinations, when newer ones from the same group are available. We compare all measurements in Figs. <ref>, <ref>, <ref>, <ref>, <ref>.

The agreement with <cit.> is in general very good, if not excellent, with the values being within 1 or 2-$\sigma$ at most, in each component (here we consider the largest of the two error-bars, since the methodology is very similar and the systematics should be directly comparable). There are however an handful of cases for which the $\mu_{\delta}$ component differs by more 2-$\sigma$: Antlia II (3.3-$\sigma$), Reticulum III (2.5-$\sigma$), Carina III (2.3-$\sigma$), Segue 1  (2.5-$\sigma$) and in principle also for some of the brightest galaxies like Sextans when only the statistical error is used. Inspection of the spatial, CM and PM location of the probable members from our code does not reveal hints of specific issues with these galaxies; the differences are likely to be the results of the methodology applied, which for some systems turns out to have a more noticeable effect.

Comparison of our systemic PM measurements (labelled "This work", shown as a black star) with literature measurements. The Gaia measurements are from Gaia Collaboration et al., 2018, Simon, 2018, Simon et al., 2020, Fritz et al., 2018, Fritz et al., 2019, Carlin & Sand, 2018, Massari & Helmi, 2018, Kallivayalil et al., 2018, Pace & Li, 2019, Pace et al., 2020, Carlin & Sand, 2018, McConnachie & Venn, 2020, McConnachie & Venn, 2020, Longeard et al., 2018, Longeard et al., 2020, Torrealba et al., 2019, Mau et al., 2020, Cerny et al., 2020, Chakrabarti et al., 2019, Gregory et al., 2020, Mutlu-Pakdil et al., 2019, Jenkins et al., 2020, Vasiliev & Baumgardt, 2021, Martínez-García et al., 2021, Li et al., 2021. Triangles indicate works that used only stars with additional information on membership, usually from  spectroscopy, but also RRLyrae stars in some cases, as Simon, 2018. HST measurements are from Piatek et al., 2003, Piatek et al., 2005, Piatek et al., 2006, Piatek et al., 2007, Pryor et al., 2015, Piatek et al., 2016, Sohn et al., 2013, Sohn et al., 2017. Among them, those indicated by diamonds (pentagons) use background galaxies (QSOs) as references. The smaller error bars include only the random Gaia error, the larger one also the systematic error when they are given as separated in the source. Note that we do not display the correlation between the PMs components here. The ellipses (when in the field of view of the plots) indicate the $\sigma=$100 km s$^{-1}$  prior (green) of McConnachie & Venn, 2020, McConnachie & Venn, 2020, and the escape velocity (grey) in the $1.6\times10^{12}$ M$_\odot$ Milky Way of Fritz et al., 2018  centred on the expected reflex motion for the system. For galaxies at a distance $>$ 500 kpc, we only plot an ellipse (purple) corresponding to 200 km s$^{-1}$.

We note that there are cases in which the statistical  MV20b uncertainties are much smaller than those we determined, i.e. in each component separately they are between 10-30% of ours for Bootes IV, Leo T, Cetus III, and as small as 2-3% for Indus I, Virgo I, Pegasus III (there are also cases in which our uncertainties are smaller, but with a reduction of at most 60%). Inspection of the output of our code corroborates the expectation of the large errors we find in these cases, given the small number of probable members with $P>0.5$ and their in general faint magnitudes; e.g. no stars with a probability of membership larger than 0.5 is found for Virgo I, Pegasus III, and only 2 for Cetus III (see "Notes on individual galaxies" for more details). MV20b find only one likely member in Cetus III, Pegasus III and Virgo I. Clearly the results for these systems need to be taken with a pinch of salt.
The inclusion or not of spectroscopy in the determination does not seem to be the culprit of this difference, since in most cases we do not notice any significant reduction in the random error when we include the l.o.s. velocity likelihood term.

On the other hand, we suspect that the main reason for the difference might be the prior in systemic PM corresponding to a 100 km s$^{-1}$ velocity dispersion used by MV20b. This can be easily seen for the three most distant galaxies of their sample (Phoenix, Eridanus II and Leo T), where they give also the motions without this prior. For the others, if we model the ratio of ours and MV20b errors as the quadratic sum of one and the ratio of our statistical error over the 100 km s$^{-1}$ dispersion, then in median our error would be only 3% larger,
with excursion in both direction our error is between 53% and 149%
of their scaled error, when we exclude Indus I (324%) and Virgo I (201%), for which we do not obtain a reliable measurement. This prior seems to be also mainly responsible for some PM differences, in cases where the absolute value of our PM is larger than in MV20b, like for Eridanus III and Horologium II.

Unlike MV20b, we do find a motion for Indus II, with 6 stars with $P>0.5$ (1 with $P>0.95$); nonetheless, the distribution of the probable members on the various planes does not transmit confidence in the result.

Typically our results compare well with the other works based on eGDR3. Li et al., 2021 start with stars that have also spectroscopic observations and then add stars without spectroscopy but with astrometric properties compatible within 5-$\sigma$ to those of the spectroscopic sample.  Their errors (which includes systematics) are nearly always larger than our total errors (in median, about 1.18$\times$,
clearly larger for many of the faintest systems). Due to their larger errors, the system with the most noticeable deviation is Bootes I (2.3$\sigma$ in R.A. when we do not apply our correction for systematics).

Martínez-García et al., 2021 focus on 14 galaxies and also use a probabilistic approach, with quality cuts likely more conservative than those adopted here. Our statistical errors are in general smaller, about 0.7x their errors; the largest differences are found for the fainter and more diffuse galaxies, probably indicating that their method needs more stars to perform well. Their systematic errors are lower, since they also model small scales effects. They also apply a QSO-based zero-point correction on their systemic PMs. Our zero-point correction and that of Martínez-García et al., 2021 differ by 0.009 mas yr$^{-1}$ on average, with the standard deviation in the ($\mu_{\alpha,*}$, $\mu_{\delta,*}$) zero-points for the galaxies in common being (0.012, 0.012) [mas yr$^{-1}$] for Martínez-García et al., 2021 and (0.013, 0.008) [mas yr$^{-1}$] in our work. When comparing our motions before QSO correction, they agree usually within 2.0 $\sigma$, the exception is Reticulum II which deviates by -2.9 $\sigma$
in R.A. (and 2.0 in Dec.) from ours and is also different from the other EDR3 determinations.

Finally,  Jenkins et al., 2020 focus on Leo IV, Leo V and Boötes I and use use only stars with spectroscopy, hence it is not surprising that their uncertainties on the systemic PMs are larger than ours. For Leo IV and Leo V the motions agree within 1 $\sigma$ and for Boötes I within 2$\sigma$ in both dimensions.

A comparison with Gaia DR2 measurements tests mainly the performance of those, due to their larger astrometric errors, but that is still useful to perform. When comparing with Fritz et al., 2018 for 38 systems in common, the standard deviation of the distribution of differences in systemic PM normalised by the uncertainty[We use only the Fritz] et al., 2018 error for normalisation, since the data from the two releases are not independent. is 0.96/1.05 in the R.A. and Dec.
component, respectively, thus within expectations. At a closer look, it appears that the standard deviation for most of the sample would be smaller than $\sim$1 and it is inflated by a few cases with larger deviations (Segue 2, 3.8; Triangulum II, 2.5 and Ursa Major I, 2.8 from the values in this work).
The accumulation at small deviation is probably understandable, since we might not be taking the correlation between the data sets well into account.

We find good agreement with the preferred values by  Fritz et al., 2019, but their sample of 4 UFDs is too small to tackle whether the small deviations found could be the result of chance.

Regarding other estimates, we note that the sub-sample of those by Kallivayalil et al., 2018 for which spectroscopy was not used do not match our motion well. That is not necessarily surprising, since in those cases Kallivayalil et al., 2018 values were based on the assumption that the galaxy was a former satellites of the LMC, which is not the case for most, like Columba I. It is more surprising for Phoenix II which likely is a LMC satellite <cit.>. Probably the reason in this case is that the LMC model used in Kallivayalil et al., 2018 does not match reality sufficiently well, as it is likely not massive enough.

The agreement is slightly worse with Massari & Helmi, 2018 with a standard deviation between ours and their results of 1.35 in R.A. and 1.06 in Dec.
This is mainly driven by Sagittarius II, which deviates by 2.9/2.0 $\sigma$ from our value, despite the large errors; nonetheless, also other of Massari & Helmi, 2018 determinations show deviations from other DR2 measurements.

As for the galaxies found beyond the MW virial radius,
we can compare with the GDR2 measurement of McConnachie et al., 2021 for NGC 6822, WLM, Leo A and IC 1613. Our statistical errors are between 27% and 47% of theirs, a clear improvement. The standard deviation between the values is 0.67/1.20
in R.A./Dec. respectively in the expected range.

Also for M33 our measurements agree reasonably well with the VLBI OH maser and HST ones by other two of Brunthaler et al., 2005, van der Marel et al., 2019.

It is interesting to compare Gaia measurements to high precision PMs obtained in a completely independent way with other telescopes like HST and VLBI. The references used by HST are either QSOs or background galaxies. The accuracy of QSOs based measurements can suffer because of the small number of reference sources, since e.g. systematic errors cannot be well derived from the data. When compared with our errors, it seems that the uncertainties quoted in the literature are underestimated, as four out of five measurements have a deviation of at least 1.8 $\sigma$ (up to 3.5 $\sigma$) in one dimension.[In contrast to the previous comparisons with Gaia measurement we apply here our QSO-based shifts, since these independent measurement are not affected by the Gaia systematics.]

In the cases with galaxies used as reference sources, 3 out of 4 works obtain deviations smaller than 1.6 $\sigma$ in both dimensions. The only exception is Sculptor (2.5-$\sigma$).
Our QSO-based shift improves the comparison for Sculptor slightly, although less than the correction adopted by Martínez-García et al., 2021; nonetheless, also with this shift a difference remains with respect to the HST determination. Since Gaia DR2 and EDR3 estimates agree with each other it seems unlikely that Gaia systematics are the only reason for it.  Nevertheless overall HST PMs based on galaxies and Gaia agree well, increasing the confidence in the precision and accuracy of both, see also the example of M31 [Salomon et al., 2020, van der Marel et al., 2012].

§ ORBIT INTEGRATION

§.§ Method

Using the PMs derived above with the distance modulus and the l.o.s. velocities from the literature listed in Tab. <ref>, we integrated the orbits of each galaxy for which spectroscopic measurements are available in three MW potentials: in two of them, the MW is treated as an isolated system (hereafter "isolated" potentials), and we explore a mass for the MW DM halo that brackets the range of likely MW masses [Boylan-Kolchin et al., 2013, Gibbons et al., 2014, Fritz et al., 2020, Wang et al., 2020]; in the other potential (hereafter "perturbed" potential) a $8.8 \times 10^{11}$  M$_{\odot}$ MW is perturbed by a $1.5 \times 10^{11}$ M$_{\odot}$ LMC, as  published by Vasiliev et al., 2021. The reason for including the latter case is that, although the mass of the LMC system is still subject to debate <cit.>, recent observations, such as the rotational velocity of the LMC [van der Marel & Kallivayalil, 2014], some perturbation of the MW's disk [Laporte et al., 2018], the dynamic of the ATLAS, Tucana III, Orphan and Sagittarius streams [Erkal et al., 2018, Erkal et al., 2019, Vasiliev et al., 2021, Li et al., 2021] and the dynamics of  distant halo stars [Erkal & Belokurov, 2020], are consistent with the idea of a massive LMC, i.e. with a mass of $1-2.5 \times 10^{11}$ M$_\odot$, perturbing significantly the gravitational potential of the MW <cit.>; therefore, we wish to investigate the impact of a massive LMC on the past orbits of the dwarf galaxies of the MW and on the account of its possible satellites.

The first isolated potential ("Light MW") is that published by Vasiliev et al., 2021 and composed of a spherical bulge of $1.2 \times 10^{10}$ M$_\odot$, an exponential disc of $5 \times 10^{10}$ M$_\odot$ and of a triaxial DM halo, with a total mass M$(<$R$_{vir})=8.8 \times 10^{11}$ M$_\odot$ within the virial radius of $251$ kpc, but unperturbed by the passage of LMC. The second isolated potential ("Heavy MW") is similar to the massive potential used by Fritz et al., 2018 and consists of a MWPotential14 [Bovy, 2015] with a more massive DM halo so that the system has a total mass of mass M$(<$r$_{\rm vir})=1.6 \times 10^{12}$ M$_\odot$ within the virial radius $r_{\rm vir} =307$ kpc .
Orbits are integrated 6 Gyrs backward and forward, with a time step of 3 Myrs, through the Agama package [Vasiliev, 2018]. The reason to integrate forward is that some of the most distant galaxies have not yet passed by pericenter, and therefore we have to integrate their orbit in the future to measure their orbital parameters.

In the perturbed potential, published by Vasiliev et al., 2021, the MW experiences the passage of a $M_{LMC}=1.5 \times 10^{11}$ M$_\odot$ LMC. In this model, the initial MW potential is as the isolated "Light MW" mentioned above, and the initial LMC is represented by a NFW halo with a scale radius of $r_s=10.84$ kpc and truncated at r$_{truc}=108.4$ kpc. In this case we integrate the orbits only backward; consider that before 5 Gyrs ago the model of Vasiliev et al., 2021 is stationary.

In order to take into account the uncertainties on the different measurements, the systematic errors in the PMs, and the correlations between the PMs along the Right Ascension and the Declination, we integrated the orbits from 100 Monte Carlo realisations of the current position and velocity of the galaxies. For the galaxies having distance modulus uncertainties lower than 0.1 mag in Table <ref>, we added in quadrature 0.1 mag, corresponding to typical systematic uncertainties between different methods of distances determinations. Additionally, we did not integrate the galaxies with heliocentric distances larger than 500 kpc since they are clearly not members of the MW system, and will require taking into account the LG in its globality <cit.>.

The current (right-handed) Cartesian Galactoncentric coordinates were calculated with Astropy [Astropy Collaboration et al., 2013, Astropy Collaboration et al., 2018], assuming that at the solar radius of $R_\odot=8.129$ kpc [Gravity Collaboration et al., 2018], the circular velocity is of $V_{circ}(R_\odot)=229.0$ km s$^{-1}$ [Eilers et al., 2019], and with a Solar peculiar motion of (U$_\odot$, V$_\odot$, W$_\odot$)=(11.1, 12.24, 7.25) km.s$^{-1}$ [Schönrich et al., 2010].

The orbital parameters, including their uncertainties, for the two isolated MW potentials are listed in Tab. <ref>. The quantities listed in this table are the Galactocentric distance of the pericenter (Peri), of the apocenter (Apo), the eccentricity (ecc), the orbit period (T), the time since the last pericenter (T$_{last,peri}$), and the fraction of galaxies reaching their apocenter in the last/next 6 Gyr $\mathcal{F}_{apo}$. For each of these quantities (except $\mathcal{F}_{apo}$), the listed values correspond to the median of that parameter found for the 100 realisation of the orbits, and the uncertainties have been calculated from the 16$-th$ and 84$-th$ quantiles. When the majority of the orbits do not reach the apocenter within the time range of the integration, the uncertainties on the apocenter, eccentricity and period cannot be computed. Therefore, for those galaxies, we rather give the value of the 16$-th$ quantile.

It is important to note here that the values given in this table are not free of biases. Indeed, the observed tangential velocity is known to be a biased estimator of the real tangential velocity of a given system [Fritz et al., 2018, van der Marel & Guhathakurta, 2008], where the former is inflated by the measurement uncertainties. Therefore, this bias also reverberates on the different derived parameters, like  apocenter or the pericenter. The bias on the tangential velocity can be understood with the following idealised experiment. Let us  assume that a galaxy is moving radially towards us with a real tangential velocity of exactly 0 km s$^{-1}$. In this case, the true pericenter is at 0 kpc. In reality, the uncertainties on the  systemic PM can only increase the observed tangential motion, this last one being positive by definition, increasing de facto the observed pericenter of the galaxy.

In order to identify the galaxies least impacted by this bias, we follow the guidelines in the Appendix of Fritz et al., 2018, where this effect was estimated with "backward" Monte-Carlo simulations. According to that study, assuming a typical velocity of $\sim$100 km s$^{-1}$ for galaxies in the MW system of satellites, one can expect the observed 3D/tangential velocity to be $<$1.5 the true one for measurements uncertainties $\lesssim$70 km $^{-1}$ on these quantities. For the rest of the analysis, we therefore mainly focus on systems with uncertainties on the total Galactocentric velocity $<$70 km $^{-1}$ (those with larger errors have their name in italics in Tab. <ref>).

§.§ Results and discussion

§.§.§ Effect of a massive LMC on the orbital properties

Patel et al., 2020 examined the changes in the orbital properties of 13 UFDs and 5 classical dSphs introduced by the inclusion of the LMC (and SMC), for a range of LMC and MW masses, using GDR2 PMs. The authors found that the orbits of both classical dSphs and UFDs were noticeably affected by the inclusion of the Magellanic Clouds, in particular the LMC for the classical dSphs;  the inclusion of the SMC did not have an impact on the number of galaxies potentially classified as satellites of the LMC, but on the longevity of this association.

Although, as we will see in Sect. <ref>, the large majority of the dwarf galaxies found around the MW that we examine were not bound to the LMC, also in this work we see that the presence of a massive LMC significantly perturbs the MW potential, modifying the orbits of most of the dwarfs. Figure <ref> gives a view on how the pericenters and apocenters are impacted within the 68% confidence interval.  In  Fig. <ref>, <ref>, <ref> we compare the orbital path of the dwarf galaxies found within 500$\,$kpc from the MW during the last 3 Gyr in the presence of the LMC ("perturbed") to the case of the isolated Light MW, which is equal to the former apart from the inclusion of the LMC and the corresponding response of the MW. This is one realisation coming from the integration of the observed systemic motions and l.o.s. velocities. For obvious reasons, the difference between the isolated and perturbed orbits become more important with lookback time since we used the same current dynamical parameters as a starting point. Thus by definition, the orbital shifting is of 0 kpc at t=0 Gyr.

Interestingly, we can see that the orbits of the majority of the galaxies are affected by the presence of the LMC, regardless of their distance to the MW or to the LMC. The influence of the inclusion of a massive LMC can manifest itself in several ways, e.g. as a change in the pericentric and/or apocentric distance, which can become larger or smaller, as well as the timing of the passages, even for systems that are very distant from the LMC and away from its orbital plane.

There are systems that are barely affected or affected in a minor way, e.g. Segue 1, Hydra II, Leo I as there are systems that are very strongly impacted, as e.g.
Draco II, Sextans, Sculptor, but even NGC 6822. For example, Sculptor would have been infalling recently onto the MW in the perturbed potential rather than being compatible with having been a long term satellite. This is mainly either due to the proximity of the LMC or to the MW reflex motion to the gravitational wake caused by the LMC infall (also called collective response)  [Garavito-Camargo et al., 2020, Petersen & Peñarrubia, 2020, Petersen & Peñarrubia, 2021, Vasiliev et al., 2021].

The perturbed potential explored here assumes a specific mass for the LMC and, of course, for the MW, but as discussed above, neither values are set in stone. That said, it is clear that if the LMC is indeed massive, the impact on the orbital properties of objects within and around the halo of the MW, and the conclusions one draw from them, can be significant. This is of course true also for individual stars in the MW (outer) stellar halo, a vast number of which will soon have 6D phase-space information thanks to large spectroscopic surveys (e.g. WHT/WEAVE, VISTA/4MOST, PSF etc.).

§.§.§ Too-big-to-fail problem/Central DM halo densities

The determination of orbital parameters of MW satellite galaxies has also been used in the literature to examine aspects of the Too-big-to-fail (TBTF) problem and make considerations on the inner DM halo densities inferred. Recently, in their analysis of PHAT-II simulations,  Robles & Bullock, 2021 found that at a given present-day maximum circular velocity, sub-haloes with small pericenters are more concentrated and have experienced a higher mass loss than those with a larger pericenter. Using GDR2 pericentric distances for the MW classical dSphs, they show that the allowed ranges for the maximum circular and peak velocities are both tightened than without the GDR2 information, with both quantities becoming smaller and going in the direction of exacerbating the TBTF problem.

In comparison with the GDR2 pericentric distances in F18, used by Robles & Bullock, 2021, those we derive here for the similar isolated potentials are tighter and towards the upper range of what the GDR2 data were suggesting for Draco and Ursa Minor (Draco: "Light MW" eGDR3 $48-56$ kpc versus 31-58kpc in GDR2 and "Heavy MW" eGDR3  $34-42$ kpc versus 21-40 kpc in GDR2), which would go in the direction of slightly alleviating the issue. Interestingly, in the perturbed potential, Draco would have a 1-$\sigma$ range of pericentric distances of 81-122kpc, while Ursa Minor would have 65-79kpc, which would push upwards the estimates of both the maximum and peak circular velocity.

Fornax is another interesting object, as dynamical modelling of the kinematic properties of its stellar component <cit.> as well as considerations and modelling of its system of globular clusters <cit.> suggest its DM halo to have a density core; this has sometimes been attributed to DM being heated up by stellar feedback <cit.> and more recently to the possibility of significant mass loss due to tides <cit.>. The range of pericentric distances determined here (see Tab. <ref>) is rather similar to that of the GDR2 determinations in F18 for the isolated potentials, probably because in the case of Fornax the uncertainty in the distance measurement plays also a role <cit.>; also the values themselves are very similar. The 1-$\sigma$ range for pericenter in the perturbed potential  (66-124kpc) is in agreement with those given by the isolated potentials. These include orbital trajectories that can significantly reduce the peak circular velocity of the DM halo due to tidal mass loss and reconcile it with the kinematic properties measured at the half-light radius Genina et al., 2020, Borukhovetskaya et al., 2021.

§.§.§ Relation to the Milky Way

Figure <ref> shows the apocentric (top) and pericentric (bottom) distances obtained in the three gravitational potentials for the systems with total Galactocentric velocity error $<$70 km $^{-1}$. Since here we wish to explore the relation to the MW, we exclude the systems that we find to be likely LMC satellites (which will be discussed in detail see below, in Sect. <ref>).

Apocentric (top) and pericentric (bottom) distance for the sample of likely MW satellite galaxies (excluding the high likely long-term LMC satellites of the LMC in Tab. <ref>) with error in 3D velocity less than 70 km s$^{-1}$. The filled and open squares show the results for the "light" and "heavy" isolated MW potentials and the light blue asterisks for the perturbed potentials. The arrows indicate those cases where either the median or the 84th percentiles were undefined. When not even the 16th percentile was defined, the symbols are placed at an apocentric distance = 900kpc. We only consider cases where the galaxy has already experienced a pericentric passage.

Uncertainties in the MW gravitational potential do cause significant variations in the orbital parameters of some of the galaxies in the sample. Nonetheless, there are some considerations that we can make.

$\bullet$ The system of Milky Way satellites From the top panel, we can see that the 16th quantile of the distribution of apocenters never reaches within the MW virial radius for Leo I and NGC 6822 within the time range of the orbit integration.
On the other hand, both galaxies do seem to have experienced one passage around the MW in the past; in particular when the perturbed and the "Heavy MW" potentials are considered for NGC 6822, this galaxy could have reached within $\sim$100kpc within the 68% confidence interval, which supports the conclusions of Teyssier et al., 2012 based on a comparison of the Galactocentric distance and radial velocity of LG field galaxies with those of haloes in the Via Lactea II simulations.

If we perform no cut in the total Galactocentric velocity error, there are other galaxies which do not have an apocenter within the MW virial radius, i.e. Eridanus II, Leo T, Phoenix, Pisces II.  However, the error in the total velocity is too large to draw meaningful conclusions without correcting for biases. Thus we cannot test yet the claim of Teyssier et al., 2012 that Leo T and Phoenix are backsplash galaxies on the basis of the orbital trajectories[Note however that the radial velocities used by those authors for Phoenix is wrong, see Kacharov] et al., 2017.  McConnachie et al., 2021 use systemic PMs to understand which ones of the isolated galaxies studied in that work might have reached within 300kpc from the MW (or M31) but our uncertainties on Leo T and Phoenix systemic motions are still too large to exclude this hypothesis in this way either. Nonetheless, the presence of young stars and HI gas in faint systems like Leo T and Phoenix supports the hypothesis that they have never approached the MW before. Also, according to McConnachie et al., 2021, the possibility that Phoenix might have entered the MW virial radius is tiny and possible only if the MW mass is at the high end of the probable range.

$\bullet$ Fast-moving galaxies Fast-moving galaxies are especially useful for placing   constraints on the MW mass  [Boylan-Kolchin et al., 2013]. A clear example is that of Leo I, whose large receding radial velocity and status as bound or unbound to the MW has caused several headaches for determinations of the MW mass since a long time <cit.>. In our determination, Leo I has a slightly smaller bulk PM value and error compared to the HST one  [Sohn et al., 2013] used by Boylan-Kolchin et al., 2013. This should cause a likely decrease in the MW mass but minor, since its Leo I total velocity is dominated by the l.o.s. component. Also the impact of the LMC does not change the main conclusions on its orbital history.

As for other fast-moving galaxies, most of those for which we can measure reliable PMs are likely to have come in with the LMC (see Sect. <ref>), as expected due its high velocity orbit, and are not of interest here.

If we concentrate on objects not classified as likely LMC satellites and that are receding, to exclude those recently infalling: in the "perturbed" potential, Bootes II would fall under this classification and has a large velocity compared to the escape velocity at its distance, as in Fritz et al., 2018; CanesVenatici I would do so in the "LightMW" but the inclusion of the LMC lowers the chances to have the apocenter beyond the MW virial radius for a 0.9$\times 10^{12}$ M$_{\odot}$ massive MW and makes the values get close to those of the "Heavy MW". The combined study of CanesVenatici I, Draco II and Hercules might turn out to be useful for considerations on the MW gravitational potential, because for the latter two objects, the likelihood of apocenters well outside of the MW virial radius has the opposite behaviour as for Canes Venatici I, i.e. it increases with the inclusion of the effect of the LMC infall.

$\bullet$ Tidal disturbances Both Bootes III and Tucana III are known to be embedded in stellar streams [Drlica-Wagner et al., 2015, Carlin & Sand, 2018].
Our analysis fully confirms the expectation that these features are the results of tidal disruption, since these two systems have pericenters in all the three potentials that bring them very close to the central regions of the MW, likely within 10 kpc or less; this is fully in line with the GDR2-based results <cit.> and the eGDR3-ones by Li et al., 2021 for Tucana III.

The stellar component of both Antlia II and Crater II has peculiar properties, with an extremely low surface brightness, large half-light radius and low l.o.s. velocity dispersion when compared to other MW satellites of similar stellar mass [Torrealba et al., 2016, Caldwell et al., 2017, Torrealba et al., 2019]. It has been postulated that also these galaxies have been "sculpted" by tidal disturbances by the MW [Fattahi et al., 2018, Sanders et al., 2018, Torrealba et al., 2019]. According to the orbital parameters derived in this work, this hypothesis appears fairly robust for Crater II, confirming the GDR2-based results, while for Antlia II it appears more sound when considering the "perturbed" and "heavy MW" potentials, than in the "light MW"  one. The pericentric distances obtained for these two models are compatible with that explored by Torrealba et al., 2019 to study whether tidal effects onto a cored DM halo could explain the low surface brightness, large half-light radius and low l.o.s. velocity dispersion of Antlia II.

The spatial distribution of the high probability member stars returned by our method shows an elongation in the outer parts for Carina's stellar component (Fig. <ref>), compatible with what seen in previous studies, based on red giant branch stars observed spectroscopically [Muñoz et al., 2006] and deep wide-area photometry [Battaglia et al., 2012, McMonigal et al., 2014]. Even if there are some intervening LMC stars in the Carina's line-of-sight, it is unlikely the feature is due to that, given that these would be included in our contamination model. Given the orbital parameters that we obtain, it appears very unlikely that this might be the result of a close interaction with the MW, nor with the LMC (see e.g. Fig. <ref>). The second last  pericenter in the perturbed potential, about 7 Gyr ago, might have brought Carina as close as 37kpc at a 1$\sigma$ level; but even if that would have been sufficient to strip its stellar component, it is highly unlikely that the elongation we see today is due to that, as tidal debris are not expected to be seen anymore after 15-20 crossing times <cit.>, i.e., between 1 and 2 Gyr in this case.

Within the 68% confidence limit, Segue 1, Segue 2, Triangulum II enter what can be potentially be seen as a dangerous zone, i.e. within 10-30 kpc from the MW centre in all the 3 potentials. We note that the clearly tidally disrupted Sagittarius (not included in this analysis) has had the most recent pericenter at about 16kpc and the second last at about 25 kpc, [Vasiliev et al., 2021], showing that tidal stripping can be efficient also at those Galactocentric distances.

Depending on the potential and/or on the error-bars some galaxies might or might not have suffered strong tidal effects (e.g. CVen I, Hercules, Willman 1, Tucana V and several more).

Especially in the perturbed potential, one can see hints that, excluding the galaxies with clearer tidal effects (streams or a very diffuse stellar component), the smallest galaxies in terms of half-light radius in general have the smallest pericenters, see Fig. <ref>. Some possible hypotheses are  that these systems had a fairly compact stellar structure at birth and have survived tides better than less compact systems,  or that this is a different manifestation of the type of features that tidal stripping might imprint on the galaxies it acts upon. A deeper investigation is deferred to the future.

Pericentric distance compared to the projected semi-major axis half light radius for the sample of likely MW  satellite galaxies (excluding the highly likely long-term LMC satellites of the LMC in Tab. <ref>) with error in 3D velocity less than 70 km s$^{-1}$. The filled and open squares show the results for the "light" and "heavy" isolated MW potentials and the asterisks for perturbed potential. The arrows indicate those cases where either the median or the 84th percentiles were undefined. We only consider cases where the galaxy has already experienced a pericentric passage.

We have investigated whether the ellipticity of the stellar component could be taken as a sign of strong tidal disturbances, but we find no clear trend between ellipticity and pericentric distances, neither in the sense of small pericenters having preferentially large ellipticities nor being preferentially round.

$\bullet$ Connection to star formation history (SFH) The dwarf galaxies that inhabit the LG have long been known to exhibit a variety of SFHs <cit.>, where this holds also when focusing only on those surrounding the MW. It is natural to ask whether a connection exists between e.g. the timing of strong enhancements of the star formation activity, or on the contrary, its quenching, and important times in the orbital history of these galaxies, e.g. infall into the MW halo or pericentric passages.

The star formation histories of the vast majority of the dwarf galaxies that surround the MW halted 8-10 Gyr ago <cit.>. At that time, in a hierarchical formation framework, the DM halo of the MW was still growing rapidly; e.g. according to the formula in Wechsler et al., 2002, around 8 Gyr ago, it would have assembled about half of its mass, while about 3 Gyr ago, already 80% of it would have been in place (this excludes the infall of LMC-like systems, which are rare in the $\Lambda$CDM cosmogony). As it can be seen e.g. in Armstrong et al., 2021, the trajectories  of MW satellites in a time evolving MW potential start deviating from those in a static potential around 3-4 Gyr ago, with the differences becoming more and more noticeable as a function of lookback time, as expected (see Fig. 4 in their article). Even though, the effect appears to be minor with respect to the mass and mass distribution of the MW, it is an additional source of uncertainty.  The time-variation of the potential is even stronger if one takes into account that around 8-10 Gyr ago, the MW accreated Gaia-Enceladus <cit.> and potentially there have been other subsequent events  <cit.>.

Therefore, since the gravitational potentials considered here do not include the growth of the MW DM halo as a function of time, we concentrate on those systems that have experienced star formation activity in the past few Gyrs and might have been linked to the MW: Leo I, Fornax, Carina, NGC 6822. We refer the reader to works such as those by Fillingham et al., 2019 and Miyoshi & Chiba, 2020 for an analysis of the connection between SFH and GDR2-based orbital trajectories of MW satellites via a comparison with simulations in the former and an analytical treatment of the growth of the MW halo in the latter. Even though both Phoenix and Leo T host young stars, we exclude them since they are currently found beyond the MW virial radius and so far there is no evidence that they might have gone through a pericentric passage. As for Eridanus II, recent SFH determinations seem to exclude the presence of young/intermediate-age stars [Simon et al., 2021, Gallart et al., 2021].

The most detailed SFHs for Leo I, Fornax and Carina were derived in Ruiz-Lara et al., 2021, Rusakov et al., 2021, de Boer et al., 2014, respectively. The star formation activity of Leo I saw the last episode of enhancement about $\sim$1 Gyr ago, after which it started decreasing till coming recently to a halt. The authors note that the timing is similar to that of the pericentric passage from HST PM measurements [Sohn et al., 2013] and GDR2 ones [Fritz et al., 2018, Gaia Collaboration et al., 2018], possibly indicating that the decrease and then stop of SFH was due to ram-pressure stripping of the gaseous component. Our eGDR3 based analysis confirms the timing of the last (and only) pericentric passage  for the 3 potentials explored (see Tab. <ref> and Fig. <ref>). Compared to Fornax and Carina, Leo I is highly likely to have come much closer to the MW center.

Within the time range considered here, Rusakov et al., 2021 detected intermittent episodes of enhanced star formation activity at $\sim$0.5, 1, 2 ago in Fornax. In our analysis, the timing of the last pericenter is at a lookback time of $-2.8_{-1.2}^{+0.8}$ Gyr (the negative sign means 2.8 Gyr ago) with an orbital period $4.4_{-0.8}^{+2.0}$ Gyr in the "Light MW" and $-1.5_{-0.4}^{+0.2}$ and an orbital period $2.5_{-0.3}^{+0.4}$ Gyr in the "Heavy MW". While it cannot be excluded that the burst occurring about 2 Gyr ago was due to a pericentric passage in either of the two isolated potentials, the orbital period would exclude that the more recent ones are due to the same cause. We examined the orbits determined for the perturbed potential and in this case, the timing of the last pericenter is similar to that of the "Heavy MW" case, with a period in the past exceeding 5 Gyr, leading to the same conclusions.

As for Carina, in the "Light MW" it has not passed pericenter yet, while in the "Heavy MW" the last pericenter occurred at a lookback time of $-1.0_{-0.3}^{+1.0}$ Gyr, with an orbital period of $2.6_{-0.5}^{+1.3}$ Gyr. The timing of the last pericenter passage in the perturbed potential is similar to that of the "Heavy MW", with the previous one occurring more than 7 Gyr ago. de Boer et al., 2014 find a long period of enhanced star formation activity, $4-6$ Gyrs ago, with a strong decrease about 1 Gyr ago, followed by an increase of SFR at the youngest ages probed, $0.25-0.5$ Gyr, till shut down. Also in this case, the decrease in star formation about 1 Gyr ago is potentially compatible with gas stripping.

Overall, a correlation can be found between some of the main features in the recent SFH of Leo I, Fornax and Carina and their orbital histories (although correlation does not necessarily implies causality!); however, the same explanation seems not to be valid for the intermittent bursts of SFH seen in the past 1-2 Gyr in these galaxies, though those events are less strong than the older star forming events.

Fusco et al., 2014 determined the SFH of NGC 6822 in 6 fields at different distances from the centre, out to 4 kpc. Their shapes differ significantly from each other and there is no unambiguous signature of a drop at a common time. On the other hand, the timing of the pericentric passage has too large uncertainties to indicate anything conclusive. In any case, given the fairly large stellar mass of NGC 6822 in comparison to the other dwarfs, and the fact that it is likely to have kept in the outer regions of the MW halo, probably the expected effect of the pericentric passage on the gas content and distribution of this galaxy should be minor and more in the direction of an outside-in ram-pressure stripping and smooth reduction of the size of the region where the bulk of the star formation has occurred, rather than sharp features in the SFH.

§.§.§ LMC satellites

Here we aim to explore which galaxies were likely part of the cortege of satellites that arrived with the LMC. For this, we used the 100 Monte-Carlo realisations of the orbits obtained in the perturbed MW potential to measure the time evolution of their relative distances to the LMC. As we can see from Fig. <ref>, the scatter increases with time. This is the consequence of the (un)accuracy of the different dynamical parameters measured for the galaxies (especially on the PMs, but also the distance), leading to a broad range of possible orbits for the galaxies with the least accurate measurements.

Although, the majority of the dwarf galaxies stay far from the LMC ($>$ 60 kpc) at every time step,  for 23 of them [Aquarius II, Böotes III, Canes Venatici II, Carina, Carina II, Carina III, Fornax, Grus II, Horologium I, Horologium II, Hydra II, Hydrus I, Phoenix II, Reticulum II, Reticulum III, Sagittarius II, Sculptor, Segue 1, Segue 2, Tucana II, Tucana III, Tucana IV, Tucana V], at least one of the 100 orbits went close enough ($<$ 60 kpc, corresponding to $\sim 6 r_s$) to potentially suggest a physical association.

For those 23 galaxies, we increased the number of Monte-Carlo realisations of the orbits to 1000, so we can measure the fraction of orbits that have been linked to the LMC in the past.
The relative position and velocity with respect to the LMC at the moment of their closest approach is listed in Tab. <ref> and is shown in Fig. <ref>. In this figure, we can see that the majority of the galaxies that pass close to the LMC are actually not linked to it, since their relative velocity is too high in comparison to the escape velocity of the LMC at any moment. On the other hand, 6 galaxies are clearly related to the LMC, with relative velocity at the moment of their closest approach significantly lower than the escape velocity: Carina II, Carina III, Hydrus I, Reticulum II, Phoenix II and Horologium I[The uncertainty on the total Galactocentric velocity of these 6 galaxies is $<$70 km s$^{-1}$.].  A closer-look to the time evolution of the distance of these galaxies from the LMC is given in Fig. <ref> and will be provided as a movie[On http://research.iac.es/proyecto/GaiaDR3LocalGroup/pages/en/audiovisual-material.php].

Interestingly, none of those galaxies are currently bound to the Magellanic system. Indeed, assuming that the gravitational potential of the LMC inside the tidal radius is unchanged by the tidal stripping of the external DM halo at any time
<cit.> and that the mass of the MW is constant, one can measure the Jacobi radius (r$_J$) as a function of time $t$, such as:

\begin{equation}
\centering
\mathrm{r}_J(t)= \mathrm{D_{LMC}}(t) \left( \frac{\mathrm{M_{LMC}}(< \mathrm{r}_{J})}{3 \mathrm{M_{MW}}(<\mathrm{D_{LMC}})} \right) \, ,
\end{equation}

where $\mathrm{D_{LMC}}$ is the Galactocentric distance of the LMC, $\mathrm{M_{LMC}}$ is the mass of the LMC inside the Jacobi radius and $\mathrm{M_{MW}}$ is the mass of the MW inside $\mathrm{D_{LMC}}$. With this formula, the Jacobi radius a t= 0, 0.5, 1.0, 1.5, 2.0, 2.5 Gyr ago is of r$_J$=19, 41, 67, 88, 104 and 106 kpc respectively. One can see that at present-day these 6 galaxies are found outside of the t=0 Jacobi radius. However, Carina II, Carina III, Hydrus I, Reticulum II and Phoenix II were still bound to the LMC at the time of closest approach, as can also be gauged by the values listed in Tab. <ref>. This allows us to conclude that they are very likely related to the LMC.

For Horologium I, the time of its closest approach to the LMC is more recent than the time when it escapes its gravitational attraction. However, over the last 5 Gyr, the distance of Horologium I relative to the LMC oscillates between 35 and 55 kpc with a velocity systematically lower than the escape velocity. Thus, we conclude that Horologium I used to be a satellite of the LMC prior to its escape. This is in agreement with the conclusions drawn by Erkal & Belokurov, 2020.

Patel et al., 2020 reached similar conclusions for Carina II, Carina  III and Hydrus I, but for Reticulum II and Phoenix II they concluded that both have been recently captured by the LMC. It is likely that the reason for these different conclusions is due to the different method used. The method used in this study is relatively similar to Erkal & Belokurov, 2020 who measured the fraction of orbits energetically bound to the LMC, while in Patel et al., 2020, they measured the orbits inside the radius of equi-density between the MW and the LMC. Moreover, the better accuracy on the systemic PMs derived in our study, which impacts significantly the probability of being related to the LMC, can also partly explain the different conclusions, as it is expected that more accurate and precise systemic PMs can deliver a stronger signal of association, when this is present, as noted by Patel et al., 2020.

As for Horologium II, despite having a median velocity relative to the LMC higher than the escape velocity, 56% of its orbits are compatible with having been recently ($<$500 Myr) ejected from the LMC system. However, it seems that Horologium II had a velocity relative to the LMC very close to the escape velocity, and this at any time. Thus, with the current precision on its systemic PM[We note that the uncertainty on its total velocity is $\sim$120 km s$^{-1}$, exceeding the quality cut of 70 km s$^{-1}$.], which dominates the uncertainties on its past orbits, it is not possible to definitively conclude if Horologium II is a former satellite of the LMC or if it has been interacting with it for a long period of time ($>$2 Gyr). A visual inspection of the different possible orbits seems to favour the first idea.

On the contrary, Grus II, despite having a relative velocity at its closest approach similar to the escape velocity, has an orbit that clearly shows that the galaxy did not originate in the Magellanic system, but has just been interacting/was captured  with/by the LMC in the last 200 Myr.

Tucana IV is also potentially parented by the LMC system. Despite its orbit reconstruction suggesting that it has been captured by the LMC about 500 Myr ago, its closest approach is at $6.57_{-2.21}^{+5.64s}$ kpc[We note that for such a small distance the positive bias likely increased the number.], 200 Myr ago with a relative velocity to the LMC lower than the escape velocity at this radius. Moreover, it has to be noticed that at that distance, the orbit of Tucana IV might have been highly perturbed by the SMC, which could have boosted its kinetic energy. Since, our model does not take into account the presence of the SMC, the orbit of Tucana IV likely overestimates its past kinetic energy, especially more than 200 Myr ago. Thus it is very likely that this galaxy has always been bound to the LMC system.

Although the majority of the possible orbits of Tucana II do not present any potential link with the LMC, 19% of its orbits have a relative velocity lower than the escape velocity of the LMC for at least one time step, while the satellite was inside the tidal radius of the LMC at that time[Hereafter we refer to this kind of orbits as "linked" orbits.] However, even in that conditions, it is very unlikely that Tucana II is related to the LMC. First because those "linked" orbits have an inclination of $\sim 47 \degree$ compared to the orbit of the LMC. Secondly, those orbits suggest that Tucana II passed through the very outer region of the halo of the LMC more than 2.8 Gyr ago, and being influenced by it $\sim 4.5$ Gyr ago. Taking into account the simplistic potential model for both the MW and the LMC on which our model is based of, it is very unlikely that Tucana II is parented to the LMC.

As for the "classical" dSphs, 25% of the potential orbits of Carina do pass through the external region of the LMC halo 2.7 Gyr ago, and for $\sim 70$% of them interacted another time with the LMC $\sim 5-6$ Gyr (see online video). Moreover, contrary to Tucana II, the orbital plane of the "linked" orbits of Carina are relatively close to the orbital plane of the LMC, with a typical angular separation between the two planes of $21 \degr$. Therefore, we cannot exclude with our study that Carina was orbiting in the external region of the LMC and has been ejected from it more than 5 Gyr, as suggested by <cit.>. However, our study tends to indicate that this scenario is unlikely, and a more accurate modelling of the MW-LMC accretion event, and/or better measurement of the current properties of Carina (especially of the distance) are required to definitively conclude something on the potential link between Carina and the Magellanic system. For Fornax the other "classical" dSphs that could be linked to the LMC, the fraction of "linked" orbits is even lower than for Carina (4%) and majority of them just pass through the LMC halo 1.8 Gyr and do not show any clear common history with the LMC, at least in the last 6 Gyr, despite having an orbital plane relatively close to the LMC, with a typical separation of of $31 \degr$. Therefore, we concluded that it is improbable that Fornax was a part of the Magellanic system.

It has to be noticed here that the fraction of "linked" orbits of Fornax and Carina that we found (0.04 and 0.25 respectively) is significantly different than the value found by Erkal & Belokurov, 2020 (0.128 and 0.004). This is the consequence of the difference in the systemic PMs that we measured above, compared to the values found with Gaia DR2 that they used for their work (see Figure <ref>).

For the two others galaxies with non zero fraction of the orbits "linked" to the LMC, Hydra II and Reticulum III, these orbits are the consequence of the large uncertainties that remain on their systemic PM, which allow a very large range of potential orbits. However, even for the few "linked" orbits, they pass only once in the external halo of the LMC in the last 5 Gyr, indicating that they are not physically associated to it.

For all the other galaxies, we can unambiguously argue that they are satellite of the MW and never used to be satellites of the LMC.

Comparison with a non perturbing LMC As said before, the mass of the Magellanic system is still heavily debated; for example based on hydro-dynamical simulations Wang et al., 2019 argued that a $10^{11}$ M$_\odot$ massive LMC cannot reproduce the Magellanic stream, neither the bridge between the LMC and the SMC. Their simulations favours a LMC with a mass of $1-2 \times 10^{10}$ M$_\odot$, which does not produce strong perturbations of the MW halo <cit.>. Thus we decided to perform the analysis of which galaxies are/were linked to the LMC also by assuming different LMC masses, ranging from $1 \times 10^{10}$ M$_\odot$ to $2 \times 10^{11}$ M$_\odot$ and assuming that it does not perturb the halo of the MW.

The result is shown on the right panel of Figure <ref>. One can see that the number of linked satellites is changing drastically depending of the mass of the LMC, with zero satellites for a LMC of $1 \times 10^{10}$, 3 (or 4 if Tucana IV is included) for a $5 \times 10^{10}$. The figure also clearly shows the importance of taking into account the perturbations produced by a massive LMC, since by neglecting them will lead to a Carina II not physically associated to the LMC even for a mass equal to $1.5 \times 10^{11}$ M$_{\odot}$, while we saw in the previous section that it is unambiguously linked to the LMC if the perturbations produced by a massive LMC are taking into account.

Estimation of the mass of the LMC
We use the number of possible LMC long-term satellites in the context of the gravitational potential perturbed by the LMC infall to estimate the mass-ratio between the Magellanic system and the MW, as done by Fritz et al., 2019. This assumes that the number of LMC/MW satellites is directly related to the mass-ratio between the two objects. Thus, the ratio of satellites between the LMC and the MW ($\mathcal{R}_{sat}$) is equal to:
$\mathcal{R}_{sat}=\frac{N_{LMC}}{N_{tot}-N_{LMC}}\, ,$
where $N_{LMC}$ is the number of LMC satellites, $N_{tot}$ is the total number of galaxies satellite, whether they are MW's or the LMC's, and $N_{MW} = N_{tot}-N_{LMC}$ the number of MW satellites. Despite being only a rough estimate, the mass so derived is still useful to verify the concordance with the mass of the LMC assumed used to find the number of LMC satellites.

Given that the exact number of $N_{LMC}$ and $N_{MW}$ is still subject to debate, mostly due to the uncertainties on systemic distances and PMs, we made 2 selections for both systems, a generous and a conservative one. For the LMC, the generous sample includes the 9 potential satellites listed in Tab. <ref> as "Highly likely parented" and "Potentially parented", plus the SMC [Murai & Fujimoto, 1980, Besla et al., 2012]; the conservative sample is composed of the 6 satellites "highly parented", plus the SMC. For the MW, we consider satellites those galaxies that have the $16^{th}$ quantile of their apocenter within the virial radius of the MW in any of the 3 potentials considered in Section <ref>, plus the Sgr dSph and the galaxies not assigned to the LMC; Sagittarius II and Crater I are not taken into account, since they are likely globular clusters <cit.>.
In the generous MW sample, we also add the galaxies for which spectroscopic measurements are not available. This leads to $N_{\rm tot}= 58$ ($N_{\rm MW}= 51-48$ for the conservative and generous LMC samples, respectively).
The conservative MW sample is restricted to the galaxies with uncertainties on their total Galactocentric velocity $< 70$ km s$^{-1}$ (see Section <ref>), which obviously excludes those without spectroscopic measurements. This leads to $N_{\rm tot}= 37$ ($N_{\rm MW}= 27-30$ for the conservative and generous LMC samples, respectively).
Doing the 4 possibles combinations that allow these 4 samples, we find a ratio of satellites between the LMC and the MW ranging from 0.14 to 0.37, with a mean of $\mathcal{R}_{sat}=0.24$, consistent with the values found by <cit.> (0.2), Erkal et al., 2018 (0.13-0.19) and Fritz et al., 2019 ($0.18^{+0.09}_{-0.08}$). This translates into a mass of the LMC between $1.5-4.1 \times 10^{11}$ M$_{\odot}$ with a mean of $2.6 \times 10^{11}$ M$_{\odot}$ for a $1.1 \times 10^{12}$ M$_{\odot}$ MW [Bland-Hawthorn & Gerhard, 2016] and between $2.0 -5.5 \times 10^{11}$ M$_{\odot}$ with a mean of $3.6 \times 10^{11}$ M$_{\odot}$ for a $1.51 \times 10^{12}$ M$_{\odot}$ MW [Fritz et al., 2020]. Although with this method the mass range for the LMC is broad, we can see that in all the case, it is consistent with a value of $\sim 1.5 \times 10^{11}$ M$_{\odot}$, but rejects the possibility of the LMC having a mass of $\sim 1-5 \times 10^{10}$ M$_{\odot}$.

Assuming an extreme scenario,  with a number of LMC satellites as low as 2 (i.e. the SMC plus one other), as we found for a low-mass non-perturbing LMC, the ratio of satellites using the conservative and generous MW samples are of $\mathcal{R}_{sat}=$0.057 and 0.036 respectively. Assuming a MW with a mass of $1.1 \times 10^{12}$ M$_{\odot}$, this correspond to a mass of the LMC of $6.27 \times 10^{10}$ M$_{\odot}$ and $3.93 \times 10^{10}$ M$_{\odot}$, respectively. This shows a lower level of consistency between the results in the case of a low mass LMC, taking into account that the mass estimated here are those for the most extreme scenario considered. Indeed, if we use the higher mass for the MW found by Fritz et al., 2020, the LMC mass estimated by its number of satellites and the mass assumed to find the number of its satellites will not be consistent anymore. However, it is important to stress that this is a ballpark estimate, since the method is based on is very simplistic. For example, it relies on the assumption that the number of luminous satellites surviving till present day is proportional to the mass of host haloes, something that is not necessarily true <cit.>. What we can say here is that the observations tend to favour a massive LMC of $1.5-2.5 \times 10^{11}$ M$_{\odot}$ rather than a $1-2 \times 10^{10}$ M$_{\odot}$ LMC.

List of galaxies reaching within 60 kpc from the LMC center for at least one out of the 100 initial Monte-Carlo realisations of their orbit in the perturbed potential. The table lists the median time of the closest approach (t$_{ca}$), the median relative distance to the LMC at that time (R$_{ca}$), the median relative velocity (V$_{ca}$) the fraction of orbit linked to the LMC for at least one time step (see text) ($\mathcal{F}_{link}$), and the median time of escape for those linked orbits (T$_{esc}$). It have to be notice that this escape time refer to the most recent moment when the galaxy escape the LMC according to the criteria we defined in the text, and do not take into account galaxy that escape the LMC and have been recaptured aftermath.

1c|t$_{ca}$ [Gyr]
1c|R$_{ca}$ [kpc]
1c|V$_{ca}$ [km.s$^{-1}$]
1cT$_{esc}$ [Gyr]

6cHighly parented to the LMC

Carina II    -0.97     $12.18_{-3.33}^{+14.02}$     $174.12_{-35.65}^{+15.77}$     1.00     -0.58

Carina III    -0.18     $13.74_{-1.54}^{+2.60}$     $164.54_{-12.12}^{+11.59}$     0.99     -0.07

Horologium I    -0.12     $36.35_{-7.20}^{+7.94}$     $113.90_{-34.31}^{+33.44}$     0.84     -0.46

Hydrus I    -0.30     $13.09_{-3.11}^{+7.14}$     $146.25_{-18.54}^{+14.48}$     1.00     -0.14

Phoenix II    -0.43     $25.99_{-14.49}^{+24.83}$     $123.57_{-47.04}^{+45.76}$     0.90     -0.31

Reticulum II    -0.19     $14.69_{-0.85}^{+1.21}$     $155.57_{-9.30}^{+7.01}$     1.00     -0.08

6cPotentially parented to the LMC

Horologium II    -0.04     $38.81_{-6.89}^{+7.11}$     $170.06_{-65.89}^{+71.39}$     0.56     -0.43

Tucana IV    -0.15     $6.57_{-2.21}^{+5.64}$     $220.88_{-45.68}^{+19.99}$     0.94     -0.06

Carina    0.00     $62.27_{-5.81}^{+5.67}$     $148.31_{-8.86}^{+10.24}$     0.25     -2.77

6cRecently captured (<1 Gyr) by the LMC

Grus II    -0.34     $25.62_{-3.33}^{+4.23}$     $186.06_{-30.59}^{+22.15}$     0.57     -0.29

6cSatellites of the MW

Aquarius II    -0.29     $34.43_{-15.08}^{+21.26}$     $396.99_{-43.99}^{+36.47}$     0.00     $-$

Bootes III    -0.20     $37.07_{-0.89}^{+1.04}$     $507.41_{-6.35}^{+3.62}$     0.00     $-$

Canes Venatici II    0.00     $196.13_{-8.30}^{+8.50}$     $380.03_{-56.31}^{+59.44}$     0.00     $-$

Fornax     -0.13     $101.33_{-7.26}^{+7.11}$     $177.90_{-20.00}^{+19.05}$     0.04     -1.81

Hydra II     0.00     $137.36_{-7.51}^{+8.02}$     $160.19_{-30.45}^{+54.30}$     0.04     -2.02

Reticulum III    -0.03     $42.96_{-11.89}^{+13.76}$     $378.09_{-69.13}^{+73.07}$     0.02     -1.43

Sagittarius II    -0.35     $18.08_{-8.55}^{+16.82}$     $269.88_{-33.21}^{+31.83}$     0.00     $-$

Sculptor    -0.11     $28.88_{-3.44}^{+3.62}$     $503.05_{-7.97}^{+7.38}$     0.00     $-$

Segue 1    -0.31     $50.61_{-5.43}^{+13.29}$     $288.02_{-68.34}^{+55.05}$     0.00     $-$

Segue 2    -0.15     $31.82_{-1.86}^{+1.89}$     $423.72_{-15.11}^{+15.66}$     0.00     $-$

Tucana II    0.00     $36.40_{-2.31}^{+3.26}$     $220.16_{-8.60}^{+9.17}$     0.19     -2.84

Tucana III    -0.08     $13.47_{-1.32}^{+1.31}$     $377.66_{-3.60}^{+3.44}$     0.00     $-$

Tucana V    -0.08     $16.69_{-8.57}^{+3.44}$     $310.79_{-18.06}^{+32.45}$     0.00     $-$

Distance from the LMC -as a function of lookback time for 100 Monte-Carlo realizations of the galaxies's orbit. The horizontal solid and dashed lines, labelled r$_s$ and r$_t$, correspond to the scale and truncation radius of the initial LMC potential, of 10.84 kpc and 108.4 kpc, respectively.

Left panel: Relative distance and velocity for the galaxies for which at least one of the 100 orbit realisations in the perturbed potential pass close to the LMC ($<$60 kpc). The black line represents the escape velocity of the initial a $1.5 \times 10^{11}$ M$_\odot$ LMC represented by a NFW profile with a scale radius $r_s=10.84$ kpc, while the dashed line shows the escape velocity of the current LMC. Each colored line shows the evolution of the relative distance and velocity of the galaxies with each point marking 100 Myr of evolution. For the galaxies that are not currently at their closest approach, the triangles show their current position on this diagram. The vertical lines at the bottom of the panel show the location of the Jacobi radius at different epochs. Right panel: Same as the left panel but assuming that the LMC follows the orbit of a point mass, and does not modify the MW potential. The black lines represent the escape velocity for a LMC with a mass of respectively $1, 5, 10, 15$ and $20 \times 10^{10}$ M$_\odot$ and a scale radius that respect the observational constraint following the requirement of Vasiliev et al., 2021.

Distance from the LMC as a function of time for the last 3 Gyr. The thicker line gives the median of the 100 random realizations and the thinner lines the 16th and 84th percentiles. The first and second row from the top depict the likely long term satellites of the LMC; the third row the possible long term satellites of the LMC and the last row the recently captured satellites.

§ FURTHER APPLICATIONS

The existence of systemic PMs and catalogues of member stars with astrometric properties opens a wealth of possibilities for the study of the internal and orbital properties of LG galaxies, which goes beyond what can be addressed in one article. Below we make a (non-exhaustive) list of the  applications that could make use either of the systemic PM measurements derived or of the list of probable member stars and associated quantities provided by this work:

* It is natural to expect that in the future LG dwarf galaxies will continue being the subject of intensive spectroscopic follow-up for the acquisition of large samples of individual stars with l.o.s. velocities and stellar atmospheric parameters (and chemical abundances). The large field-of-view and multiplex power of instruments like DESI, WEAVE, 4MOST, MOONS, PFS, MSE, and the collective power of the telescopes they are/will be mounted on, makes them particularly suitable for a comprehensive study of these systems. Lists of probable members allow to enhance the success rate of such observations by reducing the amount of contamination, as well as to assign priorities.

* The 3D motions of MW satellites, independently on their nature as galaxies or stellar clusters, can be used for determinations of the MW mass, either on their own <cit.> or as a useful addition to samples of other halo tracers, since their predominantly probe the outer parts of the MW gravitational potential and their tangential motions are known with a much higher precision than similarly distant individual MW halo stars. These 3D motions can also be used to determine the velocity anisotropy of the MW system of satellites <cit.>, typically an important ingredient for mass modeling and interesting to compare to that of other objects/tracers and to the properties of satellite systems in cosmologically simulated MW-like haloes.

* The eGDR3 3D motions of MW satellite galaxies can be used to update limits on the density of the MW hot gas corona that would be required to ram-pressure strip them of their gaseous component <cit.>. Since the common assumption is that ram-pressure stripping is most effective at pericenter, it would be interesting to perform this type of analysis both taking into account the growth of the MW DM halo (and possibly of its hot gas corona), since the SFHs of the great majority of MW satellite galaxies would suggest that they had lost their gaseous component already 8-10 Gyr ago, and combining it with models including the recent infall of a massive LMC, for those dwarf galaxies showing star formation activity in the last couple of Gyrs.

* van der Marel et al., 2019 presented the first Gaia-based study of the dynamics of the M31-M33 system, resolved the PM rotation of both galaxies, and argued that, thanks for the complete view of these rotating galaxies, the Gaia DR2 PMs determinations allow for an independent assessment of possible biases of the systemic PM measurements based on small field-of-views. The authors found that the motions of M31 and M33 support the hypothesis in which M33 is on its first infall onto M31. In our determinations of M33 motion, the statistical errors are a factor 3-4 smaller in the GDR2 based value by van der Marel et al., 2019; but systematic errors still dominate. Nonetheless, these new measurements can be used to revisit the above issues, in conjunction with eGDR3 measurements of M31 motion <cit.>.

* For those systems that enjoy both a determination of the systemic PM and l.o.s. velocity, phase-space information can be used to look for associations among dwarf galaxies, and/or with globular clusters and streams. This will be the subject of a future work. This kinematic information can also be used to determine the orbital poles to further investigates the main plane of satellites around the MW, the vast polar structure of satellites (VPOS) [Pawlowski & Kroupa, 2013]. The eGDR3 determinations  will be particularly useful for the most distant satellites, for which the PMs in DR2 were often not precise enough for a good determination of membership to the VPOS [Fritz et al., 2018, Fritz et al., 2019]. Further it is interesting to investigate how many of the members of the VPOS were once satellites of the LMC.

§ SUMMARY AND CONCLUSIONS

In this work we have jointly analyzed the spatial  distribution and the distributions onto the colour-magnitude and PM planes of  individual Gaia eDR3 sources with full astrometric solutions in the direction of 74 Local Group dwarf galaxies to determine systemic PMs of these systems. The sample includes 14 galaxies outside of the virial radius of the MW, out to $\sim$1.4 Mpc.

Our method is largely based on that by McConnachie & Venn, 2020 and McConnachie et al., 2021, to which we have introduced some modifications, aimed at a more realistic treatment of the information on the colour-magnitude diagram of the dwarf galaxies.

We are able to determine systemic PMs for 72 systems when the analysis makes no use of complementary spectroscopic information, and for 73 of them when we make use of such additional data. Overall, we consider certainly reliable the measurements for 66 of them, including all those for the galaxies outside of the MW virial radius.

The output of our analysis, including the list of members/non-members, the plots showing their  distribution on the observables used for the maximum likelihood analysis and the posterior distribution function of the systemic PMs will be made available after publication of the article.

In general our results compares very well to those in the literature based on eGDR3 data; however, we notice that our measurement uncertainties are larger than those by McConnachie & Venn, 2020 for a few systems; the main explanation for this difference is probably the prior on the velocity dispersion of the MW halo used by the other work and in a few cases differences in the quality cuts applied to the eGDR3 data.

We use eGDR3 astrometry for QSOs in the line-of-sight to the dwarf galaxies to calculate the effect of Gaia systematics on the
systemic PMs and their uncertainties. These corrections, as well as uncertainties on the distance module, are taken into account for the determination of the 3D velocities used for the reconstruction of the orbital trajectories of the galaxies around the MW (for systems out to the distance of NGC 6822).

In order to tackle the effect of the MW mass onto the orbital history and parameters, we integrate the orbits in two MW static potentials, with a mass between 0.9 and 1.6$\times 10^{12}$ M$_{\odot}$.
In addition, we complement the analysis by integrating the orbits also in a MW potential perturbed by the infall of a massive LMC, for which we use the model by Vasiliev et al., 2021. In this way, also the reflex motion imprinted onto the objects found in the outskirts of the MW halo is factored in. It should be pointed out that the errors on the transverse and 3D velocities are still very large for several systems, and this can cause biases on the determined orbital parameters. In order to limit the impact of such biases, our considerations based on the results of the orbit integration analysis mainly concern galaxies with uncertainties in the observed 3D velocities $<$70 km s$^{-1}$, which we expect to be inflated with respect to the true 3D velocity of a factor $\lesssim 0.5$.

The inclusion of a massive LMC, and the response of the MW, is found to modify the orbits of the majority of the MW satellites, regardless of their distance from the MW or the LMC, in a variety of ways, e.g. increasing/decreasing the pericentric, apocentric distance and the timing of these crucial events. Significant differences are also seen, as it would be expected, in the results from the two isolated potentials.

In general though, orbit integration of the eGDR3-based systemic PMs and literature l.o.s. velocities lead to the following conclusions in the three gravitational potentials used:

* Leo I and NGC 6822 seem not to be currently bound to the MW, although both, including NGC 6822 are likely to have entered within its virial radius once in the past.

* Bootes III and Tucana III are very likely to have reached within 10 kpc from the MW center, fully confirming the expectations that the streams in which they are embedded are the result of tidal disruption.

* The orbital properties of Crater II confirm those from GDR2-based systemic PMs and are in line with those explored by models that explain its very low surface brightness, large half-light radius and low l.o.s. velocity dispersion in the context of strong tidal disturbances from the MW.

* There are hints that the half-light radius of systems that are likely to have reached within 30 kpc from the MW center is smaller than those whose orbits kept them in more external regions.

* Concerning the SFH of MW satellites hosting young stars, i.e. Leo I, Fornax and Carina, a correlation can be found between some of the main features seen in their recent SFH, such as bursts or a clear decrease in the activity, the same explanation does not fit all the main features observed.

We also carried out an analysis aimed at identifying which ones of the galaxies surrounding the MW might have been/are physically linked to the LMC. We first identify what galaxies have at least an orbit that brings them within $\sim$60 kpc from the LMC, and then look at their total velocity and position as compared to the escape velocity curve of the LMC as a function of time. In the hypothesis that the gravitational potential and time evolution of the MW + LMC system is well represented by the Vasiliev et al., 2021 model, we find 6 systems that are highly likely to have been satellites of the LMC (Carina II, Carina III, Horologium I, Hydrus I, Phoenix II, Reticulum II), 3 that might have been potentially associated as satellites (Horologium II, Tucana IV, Carina) and one that seems to have been recently captured (Grus II). On the other hand, we find it unlikely that Fornax was associated to the LMC.

Exploring some generous and conservative estimates in the assignment of satellites to the LMC or the MW, we find that the ratio between two ranges between 0.14 and 0.37.
A simple rescaling of the number of satellite galaxies with the DM halo mass would suggest the DM halo of the LMC being in the range 1.5-4.1$\times 10^{11}$ M$_{\odot}$ for a 1.1$\times 10^{12}$ M$_{\odot}$ MW mass. It should be however pointed out that the number of (and which) galaxies are classified as LMC satellites would change when applying the same methodology using a MW potential not perturbed by the presence of the LMC and allowing for a smaller LMC DM halo mass. It would be interesting to explore further which combinations of the LMC+MW gravitational potential lead to a number of LMC satellites consistent with the expectations of cosmological theories for LMC-like haloes.

The significantly more accurate and precise eGDR3 astrometry has allowed to us to expand and improve our view of the dynamical properties of galaxies in the LG and its immediate surroundings with respect to GDR2, and even though this has been only a skimming of the potential applications of these measurements, we can only look impatiently forward to the fourth data release!

The authors acknowledge financial support through the grant (AEI/FEDER, UE) AYA2017-89076-P, as well as by the Ministerio de Ciencia, Innovación y Universidades (MCIU), through the State Budget and by the Consejería de Economía, Industria, Comercio y Conocimiento of the Canary Islands Autonomous Community, through the Regional Budget. GT acknowledges support from the Agencia Estatal de Investigación of the Ministerio de Ciencia e Innovación under grant
FJC2018-037323-I. The authors are thankful to Santi Cassisi for kindly providing the simulated Basti-IAC CMD in eGDR3 passbands.
This work has made use of data from the European Space Agency (ESA) mission
Gaia (<https://www.cosmos.esa.int/gaia>), processed by the Gaia
Data Processing and Analysis Consortium (DPAC,
<https://www.cosmos.esa.int/web/gaia/dpac/consortium>). Funding for the DPAC
has been provided by national institutions, in particular the institutions
participating in the Gaia Multilateral Agreement.

[Allende Prieto et al., 2020]
Allende Prieto, C., Cooper, A. P., Dey, A., et al. 2020, Research Notes
of the American Astronomical Society, 4, 188

[Amorisco et al., 2013]
Amorisco, N. C., Agnello, A., & Evans, N. W. 2013, , 429, L89

[Armstrong et al., 2021]
Armstrong, B. M., Bekki, K., & Ludlow, A. D. 2021, , 500, 2937

[Astropy Collaboration et al., 2018]
Astropy Collaboration, Price-Whelan, A. M., Sipőcz, B. M., et al.
2018, , 156, 123

[Astropy Collaboration et al., 2013]
Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,
, 558, A33

[Battaglia et al., 2005]
Battaglia, G., Helmi, A., Morrison, H., et al. 2005, , 364, 433

[Battaglia et al., 2008]
Battaglia, G., Helmi, A., Tolstoy, E., et al. 2008, , 681, L13

[Battaglia et al., 2012]
Battaglia, G., Irwin, M., Tolstoy, E., de Boer, T., & Mateo, M.
2012, , 761, L31

[Battaglia et al., 2015]
Battaglia, G., Sollima, A., & Nipoti, C. 2015, , 454, 2401

[Battaglia et al., 2011]
Battaglia, G., Tolstoy, E., Helmi, A., et al. 2011, , 411, 1013

[Battaglia et al., 2006]
Battaglia, G., Tolstoy, E., Helmi, A., et al. 2006, , 459, 423

[Bellazzini et al., 2014]
Bellazzini, M., Beccari, G., Fraternali, F., et al. 2014, , 566,

[Bellazzini et al., 2011]
Bellazzini, M., Beccari, G., Oosterloo, T. A., et al. 2011, , 527,

[Bellazzini et al., 2002]
Bellazzini, M., Ferraro, F. R., Origlia, L., et al. 2002, , 124,

[Belokurov et al., 2018]
Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason,
A. J. 2018, , 478, 611

[Belokurov et al., 2014]
Belokurov, V., Irwin, M. J., Koposov, S. E., et al. 2014, , 441,

[Belokurov et al., 2010]
Belokurov, V., Walker, M. G., Evans, N. W., et al. 2010, , 712,

[Belokurov et al., 2007]
Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2007, , 654, 897

[Bernard et al., 2010]
Bernard, E. J., Monelli, M., Gallart, C., et al. 2010, , 712, 1259

[Besla et al., 2012]
Besla, G., Kallivayalil, N., Hernquist, L., et al. 2012, , 421,

[Bland-Hawthorn & Gerhard, 2016]
Bland-Hawthorn, J. & Gerhard, O. 2016, , 54, 529

[Boettcher et al., 2013]
Boettcher, E., Willman, B., Fadely, R., et al. 2013, , 146, 94

[Borukhovetskaya et al., 2021]
Borukhovetskaya, A., Errani, R., Navarro, J. F., Fattahi, A., &
Santos-Santos, I. 2021, arXiv e-prints, arXiv:2104.00011

[Bovy, 2015]
Bovy, J. 2015, , 216, 29

[Boylan-Kolchin et al., 2013]
Boylan-Kolchin, M., Bullock, J. S., Sohn, S. T., Besla, G., & van der
Marel, R. P. 2013, , 768, 140

[Brown et al., 2014]
Brown, T. M., Tumlinson, J., Geha, M., et al. 2014, , 796, 91

[Brunthaler et al., 2005]
Brunthaler, A., Reid, M. J., Falcke, H., Greenhill, L. J., & Henkel,
C. 2005, Science, 307, 1440

[Caldwell et al., 2017]
Caldwell, N., Walker, M. G., Mateo, M., et al. 2017, , 839, 20

[Callingham et al., 2019]
Callingham, T. M., Cautun, M., Deason, A. J., et al. 2019, , 484,

[Cantu et al., 2020]
Cantu, S. A., Pace, A. B., Marshall, J., et al. 2020, arXiv e-prints,

[Carlin et al., 2009]
Carlin, J. L., Grillmair, C. J., Muñoz, R. R., Nidever, D. L., &
Majewski, S. R. 2009, , 702, L9

[Carlin & Sand, 2018]
Carlin, J. L. & Sand, D. J. 2018a, , 865, 7

[Carlin & Sand, 2018]
Carlin, J. L. & Sand, D. J. 2018b, , 865, 7

[Carlin et al., 2017]
Carlin, J. L., Sand, D. J., Muñoz, R. R., et al. 2017, , 154,

[Casetti-Dinescu et al., 2018]
Casetti-Dinescu, D. I., Girard, T. M., & Schriefer, M. 2018, ,
473, 4064

[Cerny et al., 2020]
Cerny, W., Pace, A. B., Drlica-Wagner, A., et al. 2020, arXiv e-prints,

[Cerny et al., 2021]
Cerny, W., Pace, A. B., Drlica-Wagner, A., et al. 2021, , 910, 18

[Chakrabarti et al., 2019]
Chakrabarti, S., Chang, P., Price-Whelan, A. M., et al. 2019, ,
886, 67

[Chambers et al., 2016]
Chambers, K. C., Magnier, E. A., Metcalfe, N., et al. 2016, arXiv
e-prints, arXiv:1612.05560

[Chiti et al., 2018]
Chiti, A., Frebel, A., Ji, A. P., et al. 2018, , 857, 74

[Cicuéndez et al., 2018]
Cicuéndez, L., Battaglia, G., Irwin, M., et al. 2018, , 609,

[Clementini et al., 2012]
Clementini, G., Cignoni, M., Contreras Ramos, R., et al. 2012, ,
756, 108

[Clementini et al., 2003]
Clementini, G., Held, E. V., Baldacci, L., & Rizzi, L. 2003, ,
588, L85

[Clementini et al., 2019]
Clementini, G., Ripepi, V., Molinaro, R., et al. 2019, , 622, A60

[Collins et al., 2017]
Collins, M. L. M., Tollerud, E. J., Sand, D. J., et al. 2017, ,
467, 573

[Conn et al., 2012]
Conn, A. R., Ibata, R. A., Lewis, G. F., et al. 2012, , 758, 11

[Conn et al., 2018]
Conn, B. C., Jerjen, H., Kim, D., & Schirmer, M. 2018a,
, 852, 68

[Conn et al., 2018]
Conn, B. C., Jerjen, H., Kim, D., & Schirmer, M. 2018b,
ArXiv e-prints [[arXiv]1803.04563]

[Conroy et al., 2019]
Conroy, C., Bonaca, A., Cargile, P., et al. 2019, , 883, 107

[Cook et al., 1999]
Cook, K. H., Mateo, M., Olszewski, E. W., et al. 1999, , 111, 306

[Correnti et al., 2009]
Correnti, M., Bellazzini, M., & Ferraro, F. R. 2009, , 397, L26

[Crnojević et al., 2014]
Crnojević, D., Ferguson, A. M. N., Irwin, M. J., et al. 2014,
, 445, 3862

[Crnojević et al., 2016]
Crnojević, D., Sand, D. J., Zaritsky, D., et al. 2016, , 824,

[Cunningham et al., 2020]
Cunningham, E. C., Garavito-Camargo, N., Deason, A. J., et al. 2020,
, 898, 4

[Dall'Ora et al., 2006]
Dall'Ora, M., Clementini, G., Kinemuchi, K., et al. 2006, , 653,

[Dall'Ora et al., 2012]
Dall'Ora, M., Kinemuchi, K., Ripepi, V., et al. 2012, , 752, 42

[Darragh-Ford et al., 2020]
Darragh-Ford, E., Nadler, E. O., McLaughlin, S., & Wechsler, R. H.
2020, arXiv e-prints, arXiv:2012.00099

[de Boer et al., 2014]
de Boer, T. J. L., Tolstoy, E., Lemasle, B., et al. 2014, , 572,

[del Pino et al., 2021]
del Pino, A., Fardal, M. A., van der Marel, R. P., et al. 2021, ,
908, 244

[Di Cintio et al., 2021]
Di Cintio, A., Mostoghiu, R., Knebe, A., & Navarro, J. 2021, arXiv
e-prints, arXiv:2103.02739

[Dolphin et al., 2003]
Dolphin, A. E., Saha, A., Skillman, E. D., et al. 2003, , 126, 187

[Drlica-Wagner et al., 2015]
Drlica-Wagner, A., Bechtol, K., Rykoff, E. S., et al.
2015a, , 813, 109

[Drlica-Wagner et al., 2015]
Drlica-Wagner, A., Bechtol, K., Rykoff, E. S., et al.
2015b, , 813, 109

[Eilers et al., 2019]
Eilers, A.-C., Hogg, D. W., Rix, H.-W., & Ness, M. K. 2019, , 871,

[Erkal et al., 2019]
Erkal, D., Belokurov, V., Laporte, C. F. P., et al. 2019, , 487,

[Erkal & Belokurov, 2020]
Erkal, D. & Belokurov, V. A. 2020, , 495, 2554

[Erkal et al., 2018]
Erkal, D., Li, T. S., Koposov, S. E., et al. 2018, ArXiv e-prints

[Errani & Navarro, 2021]
Errani, R. & Navarro, J. F. 2021, [[arXiv]2011.07077]

[Evans et al., 2018]
Evans, D. W., Riello, M., De Angeli, F., et al. 2018, , 616, A4

[Fabricius et al., 2021]