© ESO 2021

1 Introduction

This paper takes advantage and highlights the improvements from Gaia Data Release 2 (DR2) to Gaia Early Data Release 3 (EDR3) in the context of astrometry, photometry, and completeness in the Magellanic Cloud sky area. A previous Gaia DR2 science-demonstration paper on dwarf galaxies Gaia Collaboration (2018) only scratched the surface of what Gaia can tell us about these objects; it only considered their basic parameters, and barely used the photometry. Here we demonstrate how much more Gaia EDR3 shows us compared to Gaia DR2, thus demonstrating the value added by this new data release. A summary of the contents and survey properties of the Gaia EDR3 release can be found in Gaia Collaboration (2021), and a general description of the Gaia mission can be found in Gaia Collaboration (2016). Specifically, as described in Gaia Collaboration (2021), we use:

• A reduction of a factor 2 in the proper motion uncertainty.

• A new transit cross-match that provides a significant improvement in crowded areas and increases completeness.

• 33 months of data significantly reduce the Gaia scanning-law effects observed in Gaia DR2 when means and medians of parallaxes and proper motions are computed

• New photometry, with reduced systematic effects, that is less affected by crowding effects in the centre of the clouds (see Fig. 9). This helps us to unveil different stellar populations in the area of the Magellanic Clouds.

In Sect. 3 we provide an analysis of the improvements since Gaia DR2 in Gaia EDR3. In Sect. 2 we define the samples we use throughout the paper. We start by selecting objects in a radius around the centre of each cloud, and then we filter the objects using parallax, proper motions, and G magnitude. The result is two clean samples, one for the Large Magellanic Cloud (LMC) and one for the Small Magellanic Cloud (SMC). They constitute the baseline for our work. By selecting objects based on their position in the (G, G_BP − G_RP) diagram, wethen further split these samples into a set of evolutionary phase subsamples that can be used as a proxy for age selection.

In Sect. 3 we compare Gaia DR2 and Gaia EDR3 using the LMC and SMC samples. We compare the parallax and proper motion fields and show that the systematics and noise are significantly reduced. We also show that the photometry has improved by comparing the excess flux.

In Sect. 4 we use the Gaia EDR3 astrometry to resolve the 3D structure of the LMC by modelling it as a disc. We determine its parameters using a Bayesian approach. We show that the Gaia EDR3 level of parallax systematics (essentially the zero-point variations), combined with the parallax uncertainties for a distant object such as the LMC, place this determination at the very limit of feasibility. We do not reach a satisfactory result, but we conclude that it might be possible with Gaia EDR3 combined with external data, and certainly with future releases, in which the systematics and uncertainties will be reduced.

In Sect. 5 we study the kinematics of the LMC in detail. We analyse the general kinematic trends and consider the velocity profiles across the disc in detail, focusing on the separation of the rotation velocities as a function of the evolutionary stage.

In Sect. 6 we study the outskirts of the two Magellanic Clouds, and we specifically focus on one of its more prominent features: the Magellanic Bridge, a structure joining the Magellanic Clouds that formed as a result of tidal forces that stripped gas and stars from the SMC towards the LMC. We show that using Gaia EDR3 data, the Bridge becomes apparent without the need of sophisticated statistical treatment, and we can determine its velocity field and study it for different stellar populations.

In Sect. 7 we study the structure and kinematics of the spiral arms of the LMC using samples of different evolutionary phases, so that we can compare its outline as it becomes visible through different types of objects. We also study the streaming motions in the arms and produce radial velocity profiles for the different evolutionary phases. In the appendices we finally compile a variety of additional material based on Gaia EDR3 data.

2 Sample selection

We describe here the samples that we used in this paper. The selection was made in three steps that we describe below. First, we applied a spatial selection (radius around a predefined centre) to generate two base samples (LMC and SMC) in order to select objects in the general direction of the two clouds. Second, for each one of these samples, we introduced an additional selection to retain objects whose proper motions are compatible with the mean motion of each cloud. This second selection ensured that most of the contamination from foreground (Milky Way) objects was removed. Finally, we defined a set of eight subsets for each cloud based on the position in the colour-magnitude diagram (CMD) with the aim to produce groups of objects in similar evolutionary phases as a proxy of ages (see the discussion in Sect. 2.3). We did not apply the correction to G magnitudes for sources with 6p solutions that was suggested in Sect. 7.2 of Gaia Collaboration (2021). The correction is small enough (around 0.01 mag) to not have relevant effects for the methods applied in this paper, and we verified that it only very marginally affects the composition of our samples (0.04% or less of the sample size).

2.1 Spatial selection

2.1.1 LMC

The base sample for the LMC was obtained using a selection with a 20° radius around a centre defined as (α, δ) = (81.28°, −69.78°) van der Marel (2001) and a limiting G magnitude of 20.5. This selection can be reproduced using the following ADQL query in the Gaia archive:

SELECT * FROM user_edr3int4.gaia_source as g WHERE 1 = CONTAINS(POINT(’ICRS’,g.ra,g.dec), CIRCLE(’ICRS’,81.28,-69.78,20)) AND g.phot_g_mean_mag<20.5 AND g.parallax IS NOT NULL

The resulting sample contains 27, 231, 400 objects. The large selection radius causes the selection to include part of the SMC, as is shown in Fig. 1. The purpose of such a large selection area was to ensure the inclusion of the outer parts of the LMC and the regions where the LMC-SMC bridge is located.

2.1.2 SMC

The base sample for the SMC was obtained using a selection with an 11° radius around a centre defined as (α, δ) = (12.80°, −73.15°) Cioni et al. (2000a) and a limiting G magnitude of 20.5. This selection can be reproduced using the following ADQL query in the Gaia archive:

SELECT * FROM user_edr3int4.gaia_source as g WHERE 1 = CONTAINS(POINT(’ICRS’,g.ra,g.dec), CIRCLE(’ICRS’,12.80,-73.15,11)) AND g.phot_g_mean_mag<20.5 AND g.parallax IS NOT NULL

The resulting sample contains 4 709 622 objects.

2.2 Proper motion selection

Starting from the base samples described above, we followed the procedure described in Gaia Collaboration (2018) to remove foreground (Milky Way) contamination of objects based on proper motion selection. For the proper motions to be relatively easy to interpret in terms of internal velocities, we defined an orthographic projection, (see Eq. (2) from Gaia Collaboration (2018) and also Sect. 3). To determine the proper motions of the LMC and SMC and build the filters that lead to the clean samples of both clouds, we then used the following procedure. First, we computed a robust estimate of the proper motions of the clouds by:

• 1.

  We retained objects with , where r_sel is 5 deg for the LMC and 1.5 deg for the SMC.

• 2.

  We minimised the foreground contamination by selecting stars with ϖ∕σ_ϖ < 5. This parallax cut excludes solutions that are not compatible with being distant enough to be part of the LMC or SMC, and therefore possible foreground contamination from Milky Way stars. This filter was kept for the final clean samples,as described below.

• 3.

  We also introduced a magnitude limit G < 19. This limit aims to remove the less precise astrometry from the estimation of proper motions, and was relaxed to build the final clean samples, as described below.

• 4.

  We then computed median values for μ_x and μ_y with the above selection (μ_x,med, μ_y,med). These values are our reference for the typical LMC and SMC proper motions in the orthographic plane. Using these values, we determined the covariance matrix of the proper motion distribution ().

• 5.

  We retained only stars with proper motions within , where μ′ = (μ_x − μ_x,med, μ_y − μ_y,med). This corresponds to a 99% confidence region. For simplicity, we did not take the covariance matrix of individual stars into account. The aim was simply to remove clear foreground objects, and we considered the given formulation just an approximation, but sufficient for this purpose.

• 6.

  We determined the median parallax of this sample, ϖ_med, and for each star in our full sample, we determined the proper motion conditional on ϖ_med being the true parallax of the star, taking the relevant uncertainties σ and correlations ρ into account. for example, .

• 7.

  We computed new μ_x, μ_y from . We used these to repeat steps 1–4 to derive a final estimate of μ_x,med, μ_y,med, and .

Using these results, we applied the following two conditions to the base samples defined in the previous section:

• 1.

  We retained only stars with proper motions within .

• 2.

  As before, we selected only stars with ϖ∕σ_ϖ < 5 to minimise any remaining foreground contamination, but now we set a fainter magnitude limit, G < 20.5.

The resulting clean sample for the LMC contains a total of 11 156 431 objects, and the sample for the SMC contains 1 728 303 objects; their distribution in the sky is depicted in Fig. 1 and the mean astrometry is presented in Table 1. The mean parallaxes of both objects are negative, while the expected values would be (Pietrzyński et al. 2019) and (Cioni et al. 2000b). This is due to the zero-point offset in the Gaia parallaxes that was discussed in Lindegren et al. (2021a); using the values in this paper, the (rough) estimates of the LMC (− 0.0242 mas) and SMC (− 0.0185 mas) zero-points are in line with a global value of − 0.020 mas, as discussed in Sect. 4.2 of Lindegren et al. (2021a).

Fig. 1 Sky density plots for the LMC (left) and SMC (right) clean samples (after spatial and proper motion selection). Top row: plots in equatorial coordinates. Bottom row: orthographic projection (as used in Sect. 3)

2.3 Evolutionary phase subsamples

The two samples obtained following the procedure outlined in the two previous sections constitute our basic selection of objects for the LMC and SMC, our clean samples for the stars of the clouds. These were used for analysis of the LMC and SMC as a whole. A selection of basic statistics and maps using these samples is presented in Appendix A.

Several cases required a definition of subsamples that were adequate for the study of different substructures of the clouds (disc, halo, etc.), however. Ideally, we would like to select these subsamples by age, but this would require either generating our own age estimates or a cross-match with external catalogues, which is beyond the scope of a Gaia EDR3 demonstration paper such as this. Instead, we used a different approach, using a selection of samples based on the CMD of the clouds. We defined cut-outs in the shape of polygonal regions in the (G, G_BP − G_RP) diagram to select the following target evolutionary phases:

• Young 1: very young main sequence (ages < 50 Myr)

• Young 2: young main sequence ( 50 < age <   400 Myr)

• Young 3: intermediate-age main-sequence population (mixed ages reaching up 1−2 Gyr)

• RGB: red giant branch

• AGB: asymptotic giant branch (including long-period variables)

• RRL: RR-Lyrae region of the diagram

• BL: blue loop (including classical Cepheids)

• RC: red clump.

The defined areas are shown in Fig. 2. There are unassigned areas in the CMD diagrams: this is on purpose because these unassigned areas are too mixed, affected by blended stars, or too contaminated by foreground (Milky Way) stars. The areas are exclusive, that is, they do not overlap.

This rather raw selection is not even corrected for reddening, but to some extent, it can be used as an age-selected proxy. Based on a simulation using a constant star formation rate, the age-metallicity relation by Harris & Zaritsky (2009), and PARSEC1.2 models, the estimated age distribution of the resulting subsamples is shown in Fig. 3. The figure shows that the resulting subsamples indeed have different age distributions that suffice for the purposes of this demonstration paper. For the sake of brevity, we refer to these subsamples as “evolutionary phases”.

2.3.1 LMC evolutionary phases

The polygons in the CMD diagram defining the LMC subsamples are as follows, and they are represented in Fig. 2 (left panel):

• Young 1: [0.18, 16.0], [-0.3, 10.0], [-1.0, 10.0], [-1.0, 16.0], [0.18, 16.0]

• Young 2: [-1.0, 16.0], [0.18, 16.0], [0.34, 18.0], [-1.0, 18.0], [-1.0, 16.0]

• Young 3: [-0.40, 20.5], [-0.6, 19.0], [-0.6, 18.0], [0.34, 18.0], [0.40, 18.9], [0.45, 19.5], [0.70, 20.5], [-0.40, 20.5]

• RGB: [0.80, 20.5], [0.90, 19.5], [1.60, 19.8], [1.60, 19.0], [1.05, 18.41], [1.30, 16.56], [1.60, 15.3], [2.40, 15.97], [1.95, 17.75], [1.85, 19.0], [2.00, 20.5], [0.80, 20.5]

• AGB: [1.6, 15.3], [1.92, 13.9], [3.5, 15.0], [3.5, 16.9], [1.6, 15.3]

• RRL: [0.45, 19.5], [0.40, 18.9], [0.90, 18.9], [0.90, 19.5], [0.45, 19.5]

• BL: [0.90, 18.25], [0.1, 15.00], [-0.30, 10.0], [2.85, 10.0], [1.30, 16.56], [1.05, 18.41], [0.90, 18.25]

• RC: [0.90, 19.5], [0.90, 18.25], [1.60, 19.0], [1.60, 19.8], [0.90, 19.5].

The number of objects per subsample is listed in Table 2. The sky distribution of the stars in the samples is shown in Fig. A.6.

Fig. 2 Areas (as defined by the polygons given in the text) of the CMD for the LMC (left) and SMC (right) evolutionary phases. The colours are not corrected for reddening for the selection.

Fig. 3 Estimated age distribution of the selected evolutionary phase. Based on a simulation using a constant star formation rate, the age-metallicity relation by Harris & Zaritsky (2009), and PARSEC1.2 models

2.3.2 SMC evolutionary phases

The polygons in the CMD diagram defining the SMC subsamples are as follows, and they are represented in Fig. 2 (right panel):

• Young 1: [-1.00, 16.50], [-1.00, 10.00], [-0.30, 10.00], [-0.15, 15.25], [ 0.00, 16.50], [-1.00, 16.50]

• Young 2: [-1.00, 18.50], [-1.00, 16.50], [0.00, 16.50], [0.24, 18.50], [-1.00, 18.50]

• Young 3: [-0.50, 20.50], [-0.65, 20.00], [-0.65, 18.50], [0.24, 18.50], [0.312, 19.10], [0.312, 20.00], [0.50, 20.50], [-0.50, 20.50]

• RGB: [0.65, 20.50], [0.80, 20.00], [0.80, 19.50], [1.60, 19.80], [1.60, 19.60], [1.00, 18.50], [1.50, 15.843], [2.00, 16.00], [1.60, 18.50], [1.60, 20.50], [0.65, 20.50]

• AGB: [1.50,15.843], [1.75, 14.516], [3.50, 15.00], [3.50, 16.471], [1.50, 15.843]

• RRL: [0.312, 20.00], [0.312, 19.10], [0.80, 19.10], [0.80, 20.00], [0.312, 20.00]

• BL: [0.40, 18.15], [-0.15, 15.25], [-0.3, 10.00], [2.60, 10.00], [1.00, 18.50], [0.80, 18.50], [0.40, 18.15]

• RC: [0.80, 19.50], [0.80, 18.50], [1.00, 18.50], [1.60, 19.60], [1.60, 19.80], [0.80, 19.50].

The number of objects per subsample is listed in Table 3. The sky distribution of the stars in the samples is shown in Fig. A.7.

3 Comparison with DR2 results

In this section we show the improvement in astrometry and photometry of sources in the Magellanic clouds in Gaia EDR3 compared toGaia DR2. The selection of sources from Gaia DR2 for the comparison was made in the same way as for our main clean samples (as described in Sect. 2).

One of the scientific demonstration papers released with Gaia DR2, Gaia Collaboration (2018) studied the LMC and SMC, in addition to the kinematics of globular clusters and dwarf galaxies around the Milky Way. Following this study, and to ensure that the quoted (and plotted) proper motions are relatively easy to interpret in terms of internal velocities, it is particularly helpful to define an orthographic projection of the usual celestial coordinates and proper motions, (1)

where α_C and δ_C are the reference centres of the respective clouds (see Sect. 2.1).

The corresponding proper motions μ_xy = (μ_x, μ_y) and uncertainties in the form of a covariance matrix can be found from μ_α*δ = (μ_α*, μ_δ), and their uncertainty covariance matrix by the conversions (2)

where (3)

We notethat at (α_C, δ_C), we have μ_x = μ_α*, μ_y = μ_δ. We use these coordinates throughout.

In Figs. 4–7 we show the parallax and proper motion fields of the area around each of the cloud centres, as shown in the filtered Gaia DR2 and Gaia EDR3 data. We use a Voronoi binning scheme (Cappellari & Copin 2003), which produces bins with approximately 1000 stars each. The bins are therefore irregularly shaped and become large far from the centre of the clouds. Each bin is coloured according to the error-weighted mean of the indicated quantity. In each case, the dark lines are density contours.

These figures show that the Gaia EDR3 data are a clear improvement to Gaia DR2 data: the sawtooth variation that was seen in parallax and proper motion is significantly reduced. The outer bins of both the LMC and SMC still show a net positive parallax, which indicates that for these bins, foreground contamination that passes our proper motion and parallax filter makes a small but non-negligible contribution.

In Figs. 6 and 7 we show the proper motions that remain when we subtracted a linear gradient from each, so we show in each case (4)

where the central values, μ_i,0, and partial derivatives ∂μ_i∕∂x and ∂μ_i ∕∂y were evaluated as a linear fit to the values within a radius of 3° around the centre. The values found using Gaia EDR3 are shown in Table 4. This allows us to show the sawtooth pattern in proper motions more clearly. The patterns are again significantly reduced in Gaia EDR3. The faint indications of a streaming motion along the bar that were pointed out in Gaia DR2 stand out much more clearly in Gaia EDR3, and we investigate them further in Sect. 7.

As explained in Gaia Collaboration (2018, Eq. (3)), we can use the simple linear gradients to estimate the inclination, orientation and angular velocity of the disc under the assumptions that this angular velocity ω is constant, which is valid for a linearly rising rotation curve, and that the average motion is purely azimuthal in a flat disc. We define the inclination i to be the anglebetween the line-of-sight direction to the cloud centre and the rotation axis of the disc, and orientation Ω is the position angle of the receding node, measured from y towards x, that is, from north towards east. Here and elsewhere, we assume that the distances to the LMC and SMC are D_LMC = 49.5 kpc (Pietrzyński et al. 2019)and D_SMC = 62.8 kpc (Cioni et al. 2000b), respectively.

The line-of-sight velocity of the disc can either be derived from these gradients, or (as we do here) assumed given the known line-of-sight velocity of the LMC (van der Marel et al. 2002, 262.2 ± 3.4 km s⁻¹) or SMC (Harris & Zaritsky 2006, 145.6 ± 0.6 km s⁻¹). The values we find for i, Ω, and ω are 34.538°, 298.121°, 4.732 mas yr⁻¹ and 78.763°, 8.955°, 0.854 mas yr⁻¹ for the LMC and SMC, respectively. This is broadly consistent with the values found for Gaia DR2. The LMC values are consistent with those found by the more detailed investigation in Sect. 5.

In Fig. 8 we use the technique of line-integral convolution (Cabral & Leedom 1993) to better illustrate the proper motion field of the Magellanic Clouds. The direction of the lines illustrates the vector field of the proper motions, while their brightness illustrates the density (more precisely, we set the alpha parameter in MATPLOTLIB to be proportional to the 1∕4 power of the star count). The ordered rotation of the LMC is very clear from this image, while the SMC is more jumbled.

Finally, to complete this section, we compare the quality of the photometry in the LMC and SMC areas. Extracting G_BP and G_RP photometry from prism spectra is challenging in the dense, central parts of the Magellanic Clouds. A simple diagnostic for the consistency of the photometry for a source is the photometric excess factor (included in the archive), which is defined as the ratio of the flux of the prism spectra (G_BP and G_RP) and the G flux. Because the two spectra overlap slightly and have a higher response in the red, this ratio typically lies in the range 1.1–1.4 for isolated point sources, with higher values for the redder sources. Figure A.3 shows that the centres of the clouds are not very red, and Fig. 9 shows that the mean excess factor increases in these centres, but with abnormally high values in Gaia DR2 (left panel) and typical values in Gaia EDR3 (right panel). In Gaia EDR3 the background estimation has changed significantly as compared to Gaia DR2 (Riello et al. 2021), while crowding is still left uncorrected for. We conclude that the photometry in Gaia DR2 was strongly affected by background issues in the central areas, and that this problem has greatly diminished in Gaia EDR3, where traces of crowding are still visible. The G flux has only changed slightly between the two releases, that is, by a few hundredths of a magnitude, while G_BP and G_RP have been revised by a few tenths of a magnitude. It is therefore a fair assumption that the improved excess factor is driven by the improvement of G_BP and G_RP photometry in Gaia EDR3.

Fig. 4 Comparison of the parallaxes (left) and proper motions in the x and y directions (middle and right, respectively) of LMC sources in Gaia DR2 (upper panels) and Gaia EDR3 (lower panels).

Fig. 5 Sameas in Fig. 4, but for the SMC.

Fig. 6 Comparison of the residual proper motion fields of the LMC after a first-order approximation of the field was subtracted to emphasise the systematic errors in Gaia DR2 (upper panels) and Gaia EDR3 (lower panels).

Fig. 7 Sameas in Fig. 6, but for the SMC.

4 Spatialstructure of the Large Magellanic Cloud

In this section we summarise our attempts to infer the spatial distribution of sources in the LMC using a simplified model without separating the various stellar populations that constitute the galaxy. This is an oversimplification (see e.g. El Youssoufi et al. 2019, for a recent summary of the complexity of the problem when the different populations are taken into account), aimed only at exemplifying the use of the Gaia astrometry for this type of studies.

Despite the significant improvement of the Gaia EDR3 astrometry with respect to Gaia DR2, systematic problems remain, as describedin Lindegren et al. (2021b) and exemplified in the spatial distribution of median parallaxes shown in Fig. 4. In order to infer the parameters of the LMC spatial distribution, we therefore modelled the observed parallaxes as affected by a zero-point offset.

We assumed, for the sake of illustrating the magnitude of these zero-point offsets, that the sources selected as candidate members of the LMC have a mean parallax of 0.02 mas, corresponding to a distance of 50 kpc from the Sun (Pietrzyński et al. 2019). The central 90% interval around the median (binned) Gaia EDR3 parallaxes shown in Fig. 4 extends from −0.075 to 0.05 mas. We can therefore estimate the range of zero-point offsets as (−0.095,0.03). This means that the zero-point offsets are of the same order of magnitude, but larger than the expected value of the mean parallax of the LMC. Variations in parallaxes around the mean value due to the spatial distribution of the LMC sources (e.g. due to its depth or inclination angle) are expected to be much smaller. In addition, these systematics occur in combination with the usual random uncertainties associated with the individual measurements that propagate to yield the catalogue parallax uncertainty of each source. In the case of the data set used here, these parallax uncertainties have a median value of 0.17 mas. Estimating the zero-point offsets therefore is a critical element of the modelling effort we describe in this section and plays a central role in the inference of the parameters of the spatial structure of the LMC.

Unfortunately, we did not succeed in our aim of inferring geometric properties of the LMC from the Gaia EDR3 astrometric measurements. We tried several degrees of model complexity and two approaches to the inference problem: Markov chain Monte Carlo posterior sampling (MCMC) (Robert & Casella 2013), and approximate Bayesian computation (ABC) (Beaumont et al. 2002; Marjoram et al. 2003), always in the context of the Bayesian approach to inference. In the MCMC posterior sampling we used the parallaxes of the individual LMC sources to compute the full likelihood, while in the ABC approach, we binned the data in a certain number of constant-size right ascension and declination bins and employed a distance metric to compare simulations and observations in order to avoid computing the full likelihood.

Both approaches used the same probabilistic generative model for the distribution of the Gaia EDR3 parallax measurements. This model assumes that the LMC sources are spatially distributed as an elliptic double -exponential disc (similarly as in Eq. (1) of Mancini et al. 2004, but with the vertical distances from the disc mid-plane modelled by a central Laplace prior) and generates as many (proper to the disc) location coordinates as there are sources in the Gaia EDR3 sample. The model applies a number of geometrical transformations (see e.g. Weinberg & Nikolaev 2001) to generate a set of true parallaxes that are unaffected by the measurement uncertainties and/or zero-point offsets. Our generative model has nine global parameters: the disc scale length R₀, the disc scale height h₀, the disc ellipticity parameter ϵ, the disc minor-axis position angle θ_ma, and the LMC line of nodes position angle θ_LON (both angles measured with respect to the west direction), the inclination angle i of the LMC plane with respect to the plane of the sky, and the spherical coordinates of the centre of the LMC disc.

To simulate observed parallaxes, we took the Gaia EDR3 parallax uncertainties (the variance error component) and the parallax zero-point offset patterns (the systematic error component) observed in the Gaia EDR3 data into account. We modelled the latter as part of the inference process by means of a linear combination of Gaussian radial basis functions (RBFs) using the observed patterns and a canonical distance to the LMC as initial guess. Finally, each parallax measurement was simulated using a Gaussian distribution centred at the sum of the true simulated parallax and the offset generated using the RBF model.

In addition to modelling the parallax zero-point offsets using the RBF parametrisation as part of the inference process, we also tried to correct individual source parallaxes using an early version of the fit proposed in Lindegren et al. (2021b) as a function of the apparent magnitude and colour. Unfortunately, the correction is not useful for our purposes. The mentioned correction (from Lindegren et al. 2021b) is obtained by a combination of information from quasars, physical stellar pairs, and LMC sources. However, it is not able to reproduce the local variations of the parallax zero-point in the LMC field because its only dependence on positions is of the form of the sinus of the ecliptic latitude, which is almost constant in the LMC area. Additionally, the correction assumes that all the LMC stars are at the same distance embedding its internal 3D structure, which is what we aimed to determine.

In what follows we describe our attempt of using the probabilistic generative model to perform the parameter inference using the MCMC algorithm. We attempted to evaluate the full likelihood for several of the populations defined in Sect. 2. The inference process was based on a hierarchical Bayesian model and an MCMC no-U-turn posterior sampler (NUTS) (Hoffman & Gelman 2014). In this approach the true parallax of individual LMC sources was used to compute the likelihood. This implies the inclusion of one additional parameter per source (its true distance). The computational demands were so high that we were forced to distribute the likelihood computations in a TensorFlow (Abadi et al. 2016) Probability (Dillon et al. 2017) framework in the Mare Nostrum supercomputer at the Barcelona Supercomputing Centre. Unfortunately, the maturity level of the TensorFlow libraries involved was not sufficient and we did not achieve the required performance accelerations. Then, our main problem was that wewere unable to scale our models to the size of the Gaia EDR3 sample using the MCMC NUTS algorithm.

Because of the scalability issues found when using the MCMC, we decided to try with a sequential Monte Carlo approximate Bayesian computation algorithm (SMC-ABC), which is further described in Jennings & Madigan (2017) and Sect. 5 of Mor et al. (2018). The theoretical basis for these algorithms can be found in Marin et al. (2011), Beaumont et al. (2008), and Sisson & Fan (2010).

The choice of the summary statistics is crucial for the performance of the SMC-ABC algorithm. For the purposes of the present work, we defined the summary statistics as the median parallax of the stars in the LMC sample, distributed in a grid of 50 × 50 bins in right ascension (from 50° to 120°) and declination (− 50° to − 80°). The stellar sample used for this inference was the combination of the following subsamples: Young 1,Young 2,Young 3, RGB, AGB, RRL, BL, and RC.

With the SMC-ACB technique we attempted to infer up to seven parameters of the structure of the LMC: the distance to the centre D_LMC, the inclination angle i, the position angle of the line of nodes θ_LON, the position angle of the disc minor axis θ_ma, the ellipticity factor ϵ, and the position in the sky of the LMC centre . To infer these structural parameters, we chose Gaussian priors centred on the standard values found in the literature; the prior in distance is the most restrictive. Furthermore, we simultaneously inferred the model parameters of the parallax zero-point variations (i.e. the coefficients of the RBF linear model described above) using 50 basis functions. Additionally, we fixed the scale height and the radial scale length of the disc at 1.6 and 0.35 kpc, respectively.

From the SMC-ABC attempt, our conclusion is that the local parallax zero -point of the LMC in Gaia EDR3 distorts most of the signal of the 3D structure of the LMC (in the astrometry), and that there is not enough information in our summary statistics to simultaneously infer the local parallax zero-point variations and the 3D structure of the LMC. However, it may be possible if the former is constrained with additional external restrictions and/or finding an optimal way to add the information from the density distribution of the stars in the LMC area.

Fig. 8 Illustration of the proper motion field of the LMC (left) and SMC (right) using line-integral convolution. We set the alpha parameter (opacity) of the coloured lines according to the density, with the densest regions being the most opaque.

Fig. 9 Photometric excess factor, i.e. the sum of fluxes in the GBP and GRP bands over the G flux. Left: for Gaia DR2. Right: for Gaia EDR3.

5 Kinematics of the Large Magellanic Cloud

In this section we use the Gaia EDR3 data to study the kinematics of the Magellanic Clouds. The analysis is focused on the LMC because it has a clear disc structure that can be meaningfully modelled and understood; the SMC has a more complex, irregular structure that would require a more extensive and deep analysis, which is beyond the scope of this demonstration paper.

In the Sect. 5.1 we describe the method and tools we used in our analysis, and in Sect. 5.2 we present an analysis of the general kinematic trends and a detailed look at the velocity profiles in the disc, focusing on the segregation of the rotation velocities as a function of the evolutionary stage.

5.1 Method and tools

Gaia Collaboration (2018) presented formulae relating the in-plane velocities of stars to their observed proper motions under the assumption that the stars all move in a flat disc¹. Here we summarise the key results and equations.

Defining:

• a = tan i cos Ω

• b = − tan i sin Ω

• (l_x, l_y) = (sin Ω, cos Ω)

• (m_x, m_y, m_z) = (−cos i cos Ω, cos i sin Ω, sin i).

Gaia Collaboration (2018) show that Cartesian coordinates can be defined in the plane of the disc ξ, η, where (5)

and derive simultaneous equations relating the velocities to μ_x, μ_y for a given disc inclination, orientation, and bulk velocity of the galaxy. The bulk velocity of the galaxy is expressed as (μ_x,0, μ_y,0, μ_z,0), where μ_x,0 and μ_y,0 are the associated proper motions in the x and y directions at the centre of the disc, and μ_z,0 = v_z,0∕D_LMC, the associated line-of-sight velocity, expressed on the same scale as the proper motions by dividing by D_LMC. We derive (6)

Furthermore, we can gain much more physical insight by converting these Cartesian coordinates , and into polar coordinates R, ϕ, v_R, and v_ϕ.

Our strategy in this paper therefore was to fit the proper motion of the filtered LMC population as a flat rotating disc with average v_R = 0 and v_ϕ = v_ϕ(R). Our model has ten parameters, some of which can be kept fixed (based on the other knowledge of the Magellanic Clouds):

• Rotational centre of the disc on sky, parametrised as (α₀, δ₀)

• Bulk velocity in the x direction, which we parametrise as μ_x,0, the associated proper motion at the centre of the disc.

• Bulk velocity in the y direction, which we parametrise as μ_y,0, the associated proper motion at the centre of the disc.

• Bulk velocity in the z direction, which we parametrise as μ_z,0 = v_z,0∕D_LMC.

• Inclination, i.

• Orientation, Ω.

• Three parameters (v₀, r₀ and α) are used to describe the rotation curve,

To analyse the data, we considered bins of 0.08° by 0.08° in x, y in the range − 8° < x < 8°, − 8° < y < 8°. For each bin with centre x_i, y_i, we derived a maximum likelihood estimate of the typical proper motion, that is, for the ith bin, μ_i = (μ_x,i, μ_y,i), and dispersion matrix (7)

by maximising (8)

where the product is over all N_i sources in our sample in the ith bin, μ_j is the quoted proper motion of the source (μ_x,j, μ_y,j), and is the covariance matrix associated with the uncertainties as derived in Sect. 3.

We estimated the uncertainty of μ_i by bootstrap resampling within each pixel. This gave us an estimate of the error covariance matrix in proper motion for the bin, . As a simple way of taking systematic errors in proper motion into account, we added a systematic uncertainty of 0.01 mas yr⁻¹ for each component of proper motion, isotropically. This is smaller than the statistical uncertainty in most bins outside the inner ~ 3°. We chose this value because it is of the same order as the spatially dependent systematic errors found by Lindegren et al. (2021b). Binning the data allowed us to make this correction for systematic uncertainty and reduced the computational difficulty of fitting the model.

The parameters μ_x,0, μ_y,0, μ_z,0, i, and Ω give a conversion between the (x_i, y_i, μ_x,i, μ_y,i) values for each pixel and the corresponding positions and velocities in the frame of the LMC, (R_i, ϕ_i, v_R,i, and v_ϕ,i) thorugh Eq. (6). We also converted the corresponding uncertainty matrix in proper motion into one for v_R,i, v_ϕ,i (for these values of μ_x,0, μ_y,0, μ_z,0, i, and Ω), which we refer to as . The statistic that we then calculate is chi-square-like, (9)

with .

We note that the statistical uncertainties on the values we quote are very small. They are ~ 0.2 μas yr⁻¹ on μ_x,0 and μ_y,0, and less than 0.5% on the derived quantities such as μ_z,0 or i. We emphasise therefore that systematic errors, particularly those due to our simple model, are the dominant uncertainty. The difference between values in Table 5 can be seen as a gauge of these systematic errors.

Our main analysis takes the centre of rotation (α₀, δ₀) as fixed at the photometric centre α_C, δ_C (Sect. 2), and μ_z,0 = 1.115 mas yr⁻¹ taking the value from spectroscopy (Sect. 3). The parameters of this model, found by minimising χ², are given in Table 5 (along with those from the other models we considered).

We also considered the case where we did not fix μ_z,0, but left it as a free parameter. We find a value of 1.179 mas yr⁻¹, corresponding to a line-of-sight velocity of 288 km s⁻¹, which is a difference of about 7% from the value known from spectroscopy. The difference in inclination and orientation is around 1.5° in each case, and the bulk motion in x and y is almost unchanged. The ability of measuring μ_z,0 from the proper motions alone comes from the perspective contraction of the LMC as it moves away from the Sun, but we cannot expect this model-dependent result to provide a more accurate measure than from a spectroscopic study.

Finally, we considered the question of the rotational centre of the LMC. The easiest way to do this within our analysis is to allow the centre of the x, y coordinate system to shift, and then recalculate the binned values x_i, y_i, μ_x,i, and μ_y,i and uncertainties in the new coordinate system (in practice, we converted the binned values into equatorial coordinates,and then converted into the new coordinate system, rather than rebinning each time). The rotational centre of the LMC has been a matter of debate, most notably with the photometric centre and the centre of rotation for the HI gas lying at different positions. The photometric centre was found to be (81°.28, −69°.78) by van der Marel (2001), who also found that the centre of the outer isopleths (corrected for viewing angle) was at (82°.25, −69°.50). The kinematic centre of the HI gas disc has been found to be (79°.40, −69°.03) by Kim et al. (1998) or (78°.13, −69°.00) by Luks & Rohlfs (1992). Using Hubble Space Telescope (HST) proper motions in the LMC, van der Marel & Kallivayalil (2014) found a rotational centre (78°.76 ± 0°.52, −69°.19 ± 0.25°) that lies close to the centre of rotation for the HI gas, but pointed out that this was inconsistent with the rotational centre derived from studies of the line-of-sight velocity distribution in carbon stars (van der Marel et al. 2002, e.g. 81°.91 ± 0°.98, −69°.87 ± 0°.41). More recently, Wan et al. (2020) used Gaia DR2 proper motions, along with SkyMapper photometry (Wolf et al. 2018) to find dynamic centres for carbon stars, RGB stars, and young stars – (80°.90 ± 0°.29, −68°.74 ± 0°.12), (81°.23 ± 0°.04, −69°.00 ± 0°.02), and (80°.98 ± 0°.07, −69°.69 ± 0°.02), respectively.

We derive a centre of (81°.01, −69°.38) when this was left as a free parameter, which is somewhat closer to the photometric centre than to the HI centre. The inner regions of the galaxy do not provide much information in the proper motion field to find the centre because, to first order, a linearly rising rotation curve produces a linearly varying proper motion field, so that the position of the centre is degenerate with the bulk velocity. The centre of the LMC does, however, have a significant non-circular motion, which is not captured by our model, and large statistical weight in our calculations (because of the high density of stars). We therefore investigated whether cutting data from the inner few degrees of the LMC changed our results. We did this by cutting data from our analysis with (taking x and y from our original coordinate system, so that the data were the same for any centre considered), and re-deriving the parameters. The results are again listed in Table 5. The rotational centre moves slightly closer to the photometric centre as we cut larger areas from the centre of our dataset, suggesting that this result is robust against some of the incompleteness of our kinematic model. We tested whether changing the centre of our cut region affects the results (e.g. cutting data centred on the rotational centre of the HI gas instead), and the differences are very small.

In Fig. 10 we show the different proposed centres of rotation on a stellar density map of the centre of the LMC. The centres derived from Gaia EDR3 are closer to those from photometric studies than to those from the rotation of H I gas or from proper motions measured by the HST. The change in centre also naturally produces a change in derived bulk velocity, inclination, and orientation of the disc. The bulk velocity changes by ~ 0.02 mas yr⁻¹, which at the distance of the LMC corresponds to a velocity difference of 5km s⁻¹. The inclination and orientation only change by about 1°. We show plots of v_ϕ and v_R for our main model, and our model with the centre left as a free parameter (considering all data), in Fig. 11. As expected, the differences are relatively minor, although the outer parts the north-south asymmetry of the v_ϕ field is clearly reduced when the centre is left as a free parameter. The strong east-west asymmetry in v_R near the centre is also reduced (but because the minimum in v_ϕ also appearsto be offset from the centre, we are cautious about giving too much weight to this fact).

5.2 Kinematics analysis

After we robustly constrained the main parameters with the simple rotation model, we built maps of the azimuthal and radial velocities and velocity dispersions for each of the stellar evolutionary phases of the LMC, as well as for a sample combining all phases. Thislatter sample is referred to as the combined sample in this section and in Sect. 7. These maps were thus derivedat fixed and constant parameters with radius, as given by the main model of Sect. 5.1 (i = 34°, Ω = 310°, α₀ = 81.28°, δ₀ = −69.78°, μ_x,0 = 1.858 mas yr⁻¹, μ_y,0 = 0.385 mas yr⁻¹, and μ_z,0 = 1.104 mas yr⁻¹).

The angular resolution of the maps can be chosen to be as high as possible. In practice, the maps were made of 400 × 400 squared pixels of 0.04° size, which is sufficiently resolved for the simple analysis of the kinematics proposed here, and it avoids more significant statistical noise inherent to higher resolutions. At the assumed distance to the LMC, it corresponds to a linear scale of 35 pc, which is equivalent to that of observations made at 0.7′′ resolution (i.e. the typical seeing at e.g. the ESO Very Large Telescope) of galaxies located at the periphery of the local 10 Mpc volume. At this resolution, the maximum number of stars per pixel is 599, 288, 265, 239, 136, 105, 54, 52, and 13 for the combined, Young 3, RC, RGB, Young 2, BL, RRL, Young 1, and AGB samples, respectively. Despite the low surface density of the AGB, we were able to infer useful quantities, and we found that on average, AGB star kinematics compare well with other evolutionary phases. The maximum likelihood of Eq. (8) then yields the tangential and radial components of the velocity and velocity dispersion for each pixel.

Fig. 10 Positions of the centre of the LMC as found by different studies (as described in the text, SkyMapper estimates are referred to as “Sky” in the figure) superimposed on an image coloured according to the total density in star counts in the inner few degrees of the LMC. The centres found in this study lie closer to the photometric centre than to the centre of H I rotation.

5.2.1 General trends

Appendix B presents the v_ϕ, v_R, σ_ϕ, and σ_R maps for the eight stellar LMC subsamples, as well as those of the combined sample. These maps are the first of their kind ever obtained for an extragalactic disc, and the first maps that cover the integrality of the stellar kinematics for a galactic disc. To keep the description short in view of such a large quantity of kinematic data for a single galaxy, we present here example maps for two evolutionary phases only. We selected an evolved phase (RC stars) and a less evolved phase (Young 2 stars), which are both assumed to trace the kinematics of older and younger stellar ages.

The v_ϕ and v_R maps of the two phases are shown in Fig. 12 and their corresponding velocity dispersion maps in Fig. 13. They all exhibit the noisy sawtooth patterns visible in the Gaia proper motion fields (Sect. 3), as well as variations occurring at larger angular scales that may likely reflect the perturbed kinematics in the spiral arms and the bar (see also Sect. 7).

Such maps present the diversity and similarity in kinematics of the various stellar evolutionary phases. For instance, the younger phase presents higher tangential motions than the older phase (e.g. 45 versus 27 km s⁻¹ at R = 1 kpc, or 88 versus 77 km s⁻¹ at R = 4 kpc, on average), which is a beautiful signature of the asymmetric drift, while both of them present lower velocities in a region that is apparently aligned with the stellar bar, with tens of pixels sometimes at negative values (e.g. down to approximately − 25 km s⁻¹). It needs to be investigated further whether these negative values reveal counter-rotation in the bar or artificial features resulting from incorrect assumptions in this perturbed region of the disc, that is, that the stars only orbit in the z = 0 mid-plane and with v_z = 0.

The radial velocity map shows similar trends, with stronger motions for the young phase than for the RC sample. Overall, the radial motion is mostly negative for R ≲ 5 kpc, indicating inwards bulk motions towards the centre of the LMC, although this picture strongly depends on the location in the disc. Alternating negative and positive velocity patterns as a function of the azimuthal position, apparently centred on the assumed photocentre at x = y = 0, are indeed visible in the bar and spiral arms. Similarly to the rotation velocity, the velocity streaming of v_R appears to be weaker for the older stars than for the less evolved stars.

The radial and tangential velocity dispersion maps are also rich in information. Globally, the radial dispersion dominates the tangential dispersion in both samples, and the difference between the components increases with radius. There is an extended pattern of large random motions aligned with the bar in both kinematic tracers, but also a dominant feature in σ_R that is perpendicular to the bar (only for the RC sample). In this inner region of the bar, σ_R is also observed to be larger where σ_ϕ is lower for RC stars, which indicates a variation in the velocity anisotropy as a function of the azimuthal position in the bar region. As for the comparison of the samples, random motions of the RC sample are always larger than those for the young stars (e.g. 105 versus 80 km s⁻¹ in the innermost pixels, or 45 versus 20 km s⁻¹ at R = 3 kpc, on average),as expected for more evolved stars that lie in a thicker disc component than younger stars.

Fig. 11 Plots of vϕ and vR for our main model (top) and for our model with the centre left as a free parameter (bottom).

5.2.2 Velocity profiles

The 36 velocity profiles are shown in Appendix B and are made available in electronic format at the CDS. The profiles are the median values of all pixels from the maps located in radial bins of 200 pc width. This angular sampling suffices to identify variations of slope and amplitude in curves in the evolutionary phases. Radial bins with fewer than 5 pixels were discarded. The associated errors were derived from bootstrap resamplings of the velocity distributions and velocity dispersion at a given radial bin, at the 16th and 84th percentiles.

Figure 14 summarises the segregation of v_ϕ as a function of the evolutionary stage. The more evolved the stellar population, the shallower the rotation curve at low radius, and the lower the amplitude; this is an expected result from the asymmetric drift. Taking Young 1 as a reference sample with the highest values, we find that on average, the amplitude of the rotation curve of Young 2 stars is smaller by 0.6 km s⁻¹ (thus similar within the errors), and that of the BL, Young 3, AGB, RGB, RC, and RRL samples by 6, 10, 13, 17, 18, and 22 km s⁻¹, respectively. The amplitude of the combined sample lags by 15 km s⁻¹, as it is indeed dominated by the more numerous evolved stars. The BL curve is always above the Young 3 curve, and the AGB curve is above the Young 3 curve as well, but only beyond R ~ 3 kpc. Younger phases tend to have flatter rotation curves than more evolved stars. Finally, the curves of younger stars show wiggles, which are likely caused by the perturbed kinematics from the bar and spiral perturbations. The effects from the sawtooth pattern in the proper motion fields are averaged when the curves are derived, and should contribute little to the observed wiggles.

Figure 14 compares our rotation curves with a v_ϕ profile of carbon stars, as obtained by Wan et al. (2020) from modelling the Gaia DR2 astrometry. The curves of the more evolved stars from our samples agree well with their curve for R < 6 kpc. Beyondthis radius, the scatter is large in the kinematics of the carbon stars, and the curves disagree. The difference is likely caused by more significant noise in Gaia DR2 astrometry than in Gaia EDR3.

Comparisons with stellar rotation curves derived from line-of-sight velocities and HST astrometry as published in van der Marel & Kallivayalil (2014) are also shown in Fig. 14. The HST rotation curve of mixed stellar populations shown as magenta squared symbols agrees well with the Gaia curves within the quoted errors and scatter, but it has three outliers (one above 80 km s⁻¹ is not shown). The rotation velocity of old stars shown as red diamonds is systematically lower than that of our curves, while those of the young stars shown as blue circles are in fair agreement with the kinematics of the less evolved population, despite the discrepant point at R = 2.2 kpc. The large difference with the line-of-sight velocities of the older stars is not understood because the orientation parameters quoted in van der Marel & Kallivayalil (2014) do not differ strongly from those adopted here.

The v_R profiles (right panel of Fig. 14) mainly show dips with minima located at R = 2.5 − 3 kpc, near the end of the bar, except for the least evolved stars. The Young 1 and Young 2 samples indeed exhibit stronger average inwards motion at lower radius (down to v_R ~−15 km s⁻¹, R ~ 1.5 kpc). The radial motion of Young 2 stars also strongly decreases beyond R = 3 kpc. Figure 14 also shows large discrepancies between the curves of the more evolved stars with the v_R profile of carbon stars derived by Wan et al. (2020). Most of their radial velocities are > 10 km s⁻¹, and show radial motions that significantly increase as a function of radius.

Appendix B also shows the variation in the slope and amplitude of the velocity dispersion profiles as a function of the evolutionary phase. For example, the youngest phase Young 1 presents almost flat profiles, with low amplitudes (< 30 km s⁻¹), whereas the random motions of more evolved stars are steeper, and with larger amplitudes in the centre (up to 100 km s⁻¹). Again choosingYoung 1 as a reference sample, we measure that on average, σ_R of the AGB, Young 2, BL, RGB, RRL, Young 3, and RC samples is larger by 5, 21, 24, 37, 40, 40, and 52 km s⁻¹, respectively. The amplitude of the combined sample is larger by 44 km s⁻¹. Similar mean differences are observed with the tangential component of the velocity dispersion.

Fig. 12 Example of velocity maps for the LMC. Left and right columns: young and evolved evolutionary phases (Young 2 and RC, respectively, see Sect. 2.3). Upper and bottom rows: vϕ and vR, respectively.The linear velocity scales shown by colour bars are the same for the two stellar evolutionary phases and were chosen to show the structure inside the velocity fields more clearly.

Fig. 13 Same as in Fig. 12, but for the velocity dispersion.

Fig. 14 Stellar velocity curves of the LMC. Left and right panels: rotation curves and radial motions, respectively. Coloured lines are for the eight evolutionary phases and the combined sample, as given by the legend. In the left panel, blue circles, red diamonds and magenta squares are the rotation velocities for the samples of youngerand older stars from line-of-sight spectroscopic and HST astrometric measurements published in van der Marel & Kallivayalil (2014). The dashed blue lines are the tangential and radial velocity profiles derived by Wan et al. (2020) from Gaia DR2 astrometry.

6 Magellanic Bridge and the outskirts of the Magellanic Clouds

One of the most prominent features in the outskirts of two interacting galaxies is the formation of a bridge between them due to tidal forces that strip gas and stars from the least to the most massive galaxy (Toomre & Toomre 1972). The relative position of the Milky Way with respect to the Magellanic Clouds places us in the privileged position of witnessing the close encounter between the LMC and the SMC, and of studying the Magellanic Bridge.

The stellar characterisation of the structure and kinematics of the Magellanic Bridge has been pursued for a long time, with simulations (e.g. Besla et al. 2012; Diaz & Bekki 2012) and observations (e.g. Irwin et al. 1985; Harris 2007; Bagheri et al. 2013; Noël et al. 2013; Carrera et al. 2017). In addition to this expected tidally induced feature, other structures such as plumes, shells or stellar streams can be found in the outskirts of the Magellanic Clouds (e.g. Deason et al. 2017; Mackey et al. 2018; Martínez-Delgado et al. 2019; Navarrete et al. 2019). In this section we show the quality of Gaia EDR3 in highlighting the Magellanic Bridge and its kinematics, and we show several equally interesting features in the outskirts of the Magellanic Clouds.

The Magellanic Bridge was first detected as an overdensity in HI gas by Hindman et al. (1963). More recently, several studies have tried to follow the connection between the LMC and SMC using samples of stars in different evolutionary phases. Because tidal forces have similar effects on stars and gas, the Bridge would be traced by both a young stellar population with a strong correlation with the HI distribution, and an old population made of stars stripped into the Bridge by the tidal interaction. This is supported by dynamical simulations (e.g. Guglielmo et al. 2014). The stellar Magellanic Bridge was first traced by a population of young stars (Irwin et al. 1985) showing in situ star formation and a strong correlation between the location of the stars and that of HI overdensities (e.g. Skowron et al. 2014). Casetti-Dinescu et al. (2012) selected young OB-type stars in a wide area between the Clouds to study the structure and kinematics of the Bridge using GALEX, 2MASS, and the Southern Proper Motion 4 catalogue. Jacyszyn-Dobrzeniecka et al. (2020a) used Cepheids from the OGLE Collection of Variable Stars to characterise the Magellanic Bridge with young stars, while Bagheri et al. (2013) and Noël et al. (2013) used RGB stars to search for an old counterpart. Spectroscopic confirmation of stripped stars at the SMC side of the Bridge was obtained by Carrera et al. (2017). Very recently, Grady et al. (2021) assembled a catalogue of red giants from Gaia DR2 from which the authors obtained photometric metallicities using machine-learning methods. Based on the metallicity structure in the Magellanic Bridge, the authors concluded that it is composed of a mixed stellar population of LMC and SMC debris.

In this section we explore the Gaia capabilities of detecting and characterising the Magellanic Bridge using the evolutionary phase samples described in Sect. 2. Because the Magellanic Bridge encompasses the region in which the MC overlap, we have to adopt a modification of the selection described in Sect. 2. This modification takes into consideration that LMC (SMC) stars may extend farther than 20 (11) degrees and overlap with each other spatially and in proper motion. Our query is identical to the one described in Sect. 2, but we queried by HEALPIX (NSIDE=8) pixels that have a separation to their centres smaller than 35 (15) degrees from the LMC (SMC) centre. The proper motion selection described in Sect. 2 was performed, but we did not adopt any separation from either of the clouds as a membership criteria. This produced a sample that allowed for overlap in space and velocity and also provided a larger sky-coverage that is useful to explore stellar structures in the outskirts of both clouds. The total number of stars and the number in each stellar phase subsample agrees well with what we reported in Sect. 2, but because we allowed stars to mix in PM and on sky, we obtained numbers that are generally larger by <1% than in Sect. 2.

To study the Bridge, we defined two populations, one representative of the young population, and the other the RC population in both clouds. The young population was defined as the combination of Young 1_LMC, Young 2_LMC, Young 1_SMC, and Young 2_SMC, that is, inner-joined with the combination of the PM-selected LMC and SMC populations. It contains 167 643 sources. The RC population is defined in the same way, but using the RC subsamples. It contains 1 806 102 sources.

In the left panel of Fig. 15 we show a density plot of the Young stellar population in the Bridge region; the connection between the two galaxies is obvious without applying any statistical technique. The morphology of the young Magellanic Bridge is represented by an arched elongated connection between the Magellanic Clouds. In the right panel of Fig. 15 we show the density plot of the RC sample in the Bridge region. The Magellanic Bridge in the RC sample is not so clear, as expected from a more evolved and kinematically hot population, although it has been traced in RC stars by Carrera et al. (2017) at the near side of the SMC using the MAGIC spectroscopic survey. In this case, it is of key importance to remove the Milky Way foreground contamination of RC stars. This exploration has been performed by Zivick et al. (2019) using Gaia DR2 data and HST proper motions. Belokurov & Erkal (2019) used astrometry and broad-band photometry from Gaia DR2, and Schmidt et al. (2020) used data from the VISTA survey of the Magellanic Clouds (Cioni et al. 2011) and Gaia DR2 to perform a kinematic study of the region around the MC and of the Bridge region.

In Fig. 16 we specifically use the proper motions included in Gaia EDR3 to study the kinematic interaction between the Magellanic Clouds. We checked the dynamical attraction of the LMC on the SMC by plotting the vector field of the proper motion of the sources.We separately show the Young 1-2 (left panel) and the RC evolutionaryphases (right panel). In contrast to the density plot (see Fig. 15), we clearly observe, using both evolutionaryphases, a coherent motion of stars from the SMC towards the LMC. For young stars the flow moves as we would expect, from the SMC to the LMC along the Bridge (depicted in the background density plot). We emphasise that the excellent quality of the Gaia EDR3 proper motions allows tracing the interaction between theMagellanic Clouds using a rather simple strategy to separate stars into different phases of evolution.

The high quality of the Gaia EDR3 proper motions allows confirming a flow of RC stars from the SMC towards the LMC. As mentioned above, the track of an old bridge between the LMC and SMC has recently been pursued using different tracers and strategies. In this demonstration paper we considered only the RC population, which is characteristic of an intermediate to old population, and we did not use a typical > 10 Gyr old population such as that of the RR Lyrae stars. Recent works have specifically targeted the RR Lyrae stars in the bridge region of the MC (e.g. Belokurov et al. 2017; Clementini et al. 2019; Jacyszyn-Dobrzeniecka et al. 2020b). Based on their selection strategy, Belokurov et al. (2017) claimed an old RR Lyrae bridge for Gaia DR1 RR Lyraes. The Gaia DR2 bona fide RR Lyraes and those from the Gaia EDR3 sample (see Sect. 2) both show a smooth halo-like density distribution (Clementini et al. 2019), however. The Gaia DR2 accompanying paper was confirmed by Jacyszyn-Dobrzeniecka et al. (2020b) using the extended OGLE catalogue. Evans et al. (2018) stated in a Gaia DR2 accompanying paper that a suboptimal computation affected the mean magnitude standard deviation given in Gaia DR1 and Gaia DR2 (and revised in Gaia DR3; Riello et al. 2021), which may have affected the selection strategy of Belokurov et al. (2017) with only candidate RR Lyrae stars. We show here that a flow of RC stars (see Fig. 16) confirms a bridge composed of an evolved population, and it would have a similar trajectory to that of Belokurov et al. (2017). It is beyond the scope of this paper, however, to make a quantitative comparison.

In addition to the Bridge, de Vaucouleurs & Freeman (1972a) showed a wealth of substructure in the outskirts of the Magellanic Clouds in the 1970s. More recently, new shells, plumes, and streams have been detected using different surveys or photometric techniques (e.g. Pieres et al. 2017; Belokurov & Erkal 2019; Martínez-Delgado et al. 2019). To search for substructures around the Clouds, we adopted a more restrictive selection using the RGB and RC subsamples. First, we corrected for foreground extinction using Schlegel et al. (1998) (with the correction from Schlafly & Finkbeiner 2011), and we adopted a (Cardelli et al. 1989) extinction curve with R_V = 3.1. This correction is accurate in the outskirts of the Clouds because there is little internal extinction from the LMC and SMC themselves. Second, we built a tighter colour-magnitude selection polygon based on the extinction-corrected RC and RGB samples; in this case, we are stricter in the colour range allowed for these two evolutionary phases as in the RC and RGB samples described in Sect. 2. Additionally, we applied a magnitude cut of G < 19 and selectedonly stars with a parallax smaller than 0.15. This led to a sample of stars that is less strongly affected by Milky Way foreground, thus allowing us to explore faint substructures in the outskirts; we call this selection LMCout.

In Fig. 17 we show a star count map of the Magellanic Cloud region to highlight the substructure found using Gaia EDR3. We also annotate a few notable features and show the measured total dispersion and velocity vector map (velocities and dispersion were computed using LMC-centred coordinates). The northern tidal arm (NTA) reported initially by Mackey et al. (2016) is visible in the figure, and this structure is also visible as a velocity low-dispersion feature, with velocities consistent with the LMC main body. A southern tidal arm (STA) (Belokurov & Erkal 2019) is also evident, which shows indications of being dynamically cold, like the NTA, and the velocities are consistent with those of the LMC. The STA appears to be connected with the SMC through a narrow elongation east of the SMC. The SMC northern overdensity (Pieres et al. 2017) is also evident, and a spatially thinner structure is also seen to be connected to it. In addition to these known substructures, we find a faint overdensity east of the LMC that is visible in the velocity field and density map. We also see a similar structure, but more conspicuous, on the western side of the LMC, close to y_LMC = 0. We note, however, that the features observed near y_LMC = 0 coincide with a region of elevated number of Gaia transits. The eastern feature is also prominent in near-infrared maps from the VISTA Hemisphere Survey (El Youssoufi et al. 2020, priv. comm.).

Fig. 15 Sky density plot for the Bridge region using the Young1 and Young2 evolutionary phases (left panel) and the RC evolutionary phase (right panel). The bin size is 0.06 deg in x and y. The colour bar is in log scale and the black contours are at the levels 0.1, 1, 5, and 70 (young evolutionary phases) and 0.25, 1.5, 15, and 165 (RC sample).

Fig. 16 Vector field of the proper motions in the Magellanic Clouds using theYoung 1 +Young 2 (left panel) and RC (right panel) samples. The coordinates are centred on the SMC. In the background, to guide the eye of the reader, we show the density in logarithmic scale. The dashed line in the left panel shows the location of the young bridge from Belokurov et al. (2017). The velocity vector field is only shown for bins with more than 10 (200) stars in theYoung 1 +Young 2 (RC) sample.

Fig. 17 Number count map for the LMCout sample (left): the names of the most notable substructures are given and, for a better visualisation, the inner parts of the LMC and SMC were omitted. Total velocity dispersion map (right): the velocity vector field is shown as orange arrows.

7 Spiral arms in the Large Magellanic Cloud

The LMC is a prototype of barred Magellanic spiral galaxies that are characterised by an off-centre bar and one prominent spiral arm. The dynamical interactions between the LMC and the SMC (e.g. Besla et al. 2016) are probably responsible for this and other spiral features associated with the galaxy (de Vaucouleurs & Freeman 1972a). A comprehensive study of the morphology of the LMC based on near-infrared observations is given in El Youssoufi et al. (2019), where high spatial resolution maps (0.03 deg²) of stellar populations with different median ages show at least four distinct spiral features. These arms emerge predominantly from the ends of the bar, one in the east extending to the south, and three in the west, one of which extends north, one north-west (the most prominent arm), and the third extends south. The arms are well traced by stellar populations younger than a few million years, while old stellar populations instead show external features that may be associated with a ring-like structure (e.g. Choi et al. 2018). The long-term stability of the prominent spiral arm was studied by Ruiz-Lara et al. (2020) using deep optical photometry to derive the star formation history throughout the galaxy. This structure could have formed a few million years ago at the time when the Magellanic Stream and the Leading Arm formed as well from a close encounter between the LMC and the SMC. The authors concluded that the distribution of HI gas and the coherent star formation at the location of the arm support this scenario. In this section we show that the spiral structure of the LMC can be highlighted and studied with the Gaia EDR3 data.

7.1 Basic properties as a function of evolutionary phase

We discuss here the appearance of the spiral arms of the LMC using the evolutionary phase samples as proxies for age-selected samples (see Sect. 2.3). In Fig. 18 we show the maps obtained in the LMC for these samples. The maps were constructed considering a region 20 × 20 deg² around the galaxy, applying a Gaussian smoothing and sampling with 400 × 400 pixels, each with a dimension of 3′ × 3′. Each map was normalised to the total number of objects. The figure shows that the main structures of the LMC, that is, the bar and the spiral arms, are clearly outlined by BL stars, objects with ages in the range of 50–350 Myr. We therefore chose this population as a reference for the comparison of the spiral structure(s) in other stellar populations of different ages.

Because the differential maps of the BL with respect to Young 1 andYoung 2 were similar, we merged these two evolutionary phases into a single sample. We refer to this merged sample as the young population of stars with age < 400 Myr, which is shown in the middle top panel of Fig. 18. Similar considerations applied to the RC, RGB, and RRL populations, all older than 1 Gyr, and we refer to these merged evolved populations as the old population. The relative differential map is shown in the lower panel of Fig. 18. Finally, the differential map with respect to the Young 3 population (MS stars with ages < 1−1.5 Gyr) is shown in the middle panels of the same figure.

The analysis of the top panel in Fig. 18 reveals that the young population is more concentrated around the bar and an inner northern arm, showing a clumpy structure. The residual map with respect to the BL shows that this latter population has an excess of stars along the bar, in the spiral feature at the end of the bar, and in an outer north-east arm (referred to as the eastern arm hereafter).

The comparison of BL andYoung 3 populations shows that the two populations are distributed in a very similar way, even though the BL still displays an excess along the bar, especially in the eastern region, where it shows a concentration superior to any other population in the LMC (ΔRA, ΔDec ~ −3°, −1°). The older populations of the LMC have a homogeneous distribution; the star density decreases smoothly from the centre to the outskirts of the LMC. The lower density along the bar is caused by the Gaia incompleteness in this crowded region. The difference with the BL population again shows an excess of stars along the bar and the above-mentioned clump in the eastern bar, but these features might in part be justified by the incompleteness of old populations in the more central region of the bar. In contrast, the excess of BL stars in the inner and outer arms appears to be genuine.

7.2 Strength and phase of the density perturbations

To be more quantitative on the effect of the bar and spiral structure on the stellar density, we performed discrete fast Fourier transforms (FFTs) of density maps of the BL and the combined samples. For this purpose, we again used 400 × 400 pixels maps, but at 0.04° sampling, as in Sect. 5.2. This allowed us to estimate the properties of any asymmetries in the density maps.

Because the apparent dominant modes of perturbations are the bar (second-order perturbation), the inner spiral structure starting at the end of the bar, and the eastern outer arm (first-order perturbation), we present results up to the second-order harmonics, although the discrete FFTs yield as many orders as existing pixels in a vector. Therefore the analytic form equivalent to the discrete FFT applied to a density map is Σ(R, ϕ) = ∑_kΣ_k(R)cosk(ϕ − ϕ_k(R)), where k is an integer, ϕ the azimuthal angle in the plane of the LMC with the reference ϕ = 0 chosen aligned to the photometric major axis of the disc (Ω_phot), Σ₀ the axisymmetric surface density, and Σ_k and ϕ_k are the amplitude and phase of the kth asymmetry.

We measured Ω_phot by isophotalellipse fitting to the stellar surface density map of the combined sample. To avoid confusion with Sect. 5, which gives radii in a kinematic frame oriented along the kinematic position angle of Ω = 310°, we refer to R_phot as the galactocentric radius measured in the photometric frame, which is aligned on Ω_phot. With this, we find the bar semi-major axis at a position angle of ~105°, and define the one for the disc as that of the average value in the radial range R_phot = 3.5−7.5 kpc, which is Ω_phot ~ 10°. The photometric major axis therefore differs by ~60° from the kinematic major axis. A similar discrepancy has been reported in van der Marel et al. (2002).

Figure 19 shows the results of the Fourier transforms. We restrict the analysis to R_phot ≤ 7.5 kpc. The axisymmetric density profile of BL stars is more centrally peaked than that of the whole sample. At the peak of Σ₁ and Σ₂, the strengths of the lopsided outer spiral arm (k = 1) and bisymmetry (k = 2) reach 60% and 40% the amplitude of the axisymmetric mode for the combined sample, and 48 and 40% for the BL stars. At R_phot > 5.5 kpc, the strengthof the lopsided spiral is similar to the axisymmetric value. The lopsidedness and bisymmetry perturbations are therefore not negligible in the LMC.

Both samples show that the dominant perturbation at low radius is the bisymmetric mode (R_phot < 3.5 kpc for BL stars, R_phot < 2.5 kpc for the combined stars), while the lopsided mode dominates at larger radii. In the inner kpc, the Σ₂ profile of the whole sample presents a dearth of stars that is lacking in the density map of younger stars. This is caused by the incompleteness of Gaia in this crowded region of the LMC disc. This dearth of stars also affects the inner profile of ϕ₂ for all stars as a central phase dip.

The orientation of the inner k = 2 perturbation does not change much in the inner disc, with a bar oriented with a phase angle of 1.6π rad (modulo π) for both samples. The k = 2 spiral structure of BL stars has a phase angle of 1.4π rad (R_phot > 3 kpc), while that for the combined sample smoothly increases to ~2.1π rad for R_phot = 3.5−5.2 kpc, then remains constant out to the last radius. The phase angles of the lopsided mode continuously vary as a function of radius, and the two stellar samples present different shapes of ϕ₁. The similar shape of ϕ₁ and ϕ₂ at R_phot > 2.5 kpc for the combined sample of stars is remarkable, and the amplitudes only differ by less than 0.2π. The outer spiral structure in the LMC combined sample is thus made of two modes that are tightly coupled.

Fig. 18 Differential density maps for different stellar populations in the LMC. Left column: density map in the BL evolutionary phase. Middle panels from top to bottom: density maps of the young phase, i.e. stars with age <400 Myr, Young 3 phase, i.e. main- sequence stars with age 1−1.5 Gyr, and oldphase, i.e. stars with age >1 Gyr. Right column: normalised difference between the left and middle panel maps, shown with a logarithmic stretch.

7.3 Across and along streaming motions in the eastern spiral arm

The LMC velocity fields have been shown to exhibit variations stemming from the juxtaposition of an observational sawtooth-like pattern inherent to Gaia, and others likely caused by intrinsic perturbations of the gravitational potential of the LMC (Sect. 5.2). Here we illustrate in a simple way the variation of v_ϕ and v_R along and across the dominant outer k = 1 spiral arm to the east in the combined sample.

To isolate the effects of the outer arm better, we only considered the region where the inner k = 2 mode becomes negligible, that is, all pixels located at R_phot > 3.8 kpc (Fig. 20). We built azimuth-radius diagrams of the stellar density and tangential and radial velocities by calculating average star counts, v_ϕ and v_R in bins of 5° size in azimuthal angle, and 63 pc in radius (Fig. 21). The horizontal variation is thus a good proxy of the streaming of v_ϕ and v_R along the eastern spiral arm, while the vertical axis is a good proxy for the velocity variation across the spiral arm.

The density of the spiral arm is strongly asymmetric as a function of azimuthal angle, caused by its lopsided nature. The uppermost isocontour of density (mean star count of ~ 20 stars) approximately delineates the maximum radial extent of the spiral arm, which extends to R_phot ~ 7 kpc along the photometric major axis (ϕ = 0°) to R_phot ~ 5.4 kpc (ϕ = 240°). The highest densitiesaround ϕ = 300° at lower radii correspond to regions of the LMC that are part of the inner spiral structure, thus not strictly belonging to the outer lopsided eastern arm.

Along the horizontal axis, v_ϕ is maximum in higher density regions and minimum in lower density regions. When we consider pixels below the outermost contour, the azimuthal streaming in the arm is relatively constant (60 ≲ v_ϕ ≲ 90 km s⁻¹). An exception to this occurs at R ~ 6 kpc owing to the lower values of v_ϕ around ϕ = 100°. As the pixels above the uppermost contour likely probe stars beyond the spiral arm, the difference in colours between pixels below (redder) or above (bluer, v_ϕ < 60 km s⁻¹) the uppermost contour clearly shows the effect of the arm on v_ϕ in the azimuthal direction. The radial velocity also varies significantly with azimuth. It is stronger in higher density regions around ϕ = 50° and ϕ = 200° (v_R > −10 km s⁻¹) and in lower density regions for 100 < ϕ < 180°, but with inward motions (v_R < −10 km s⁻¹). The noise in v_R is higher outside the arm at large radii.

Along the vertical axis, v_ϕ is observed to decline with radius across the spiral arm, and the decrease is not complete at the same rate for different azimuthal angles. This implies a wide diversity of shapes and amplitudes in the LMC rotation curve as a function of azimuth. We have observed this trend in the v_ϕ map of Sect. 5.2. The radial velocity also varies across the spiral arm, but there appears to be no clear rule, unlike for v_ϕ. For example, the peak of v_R at ϕ = 50° occurs at R_phot = 5.5 − 6 kpc, thus beyond the location of the density peak(R_phot ≲ 4.5 kpc). However, at an angle of for instance ϕ = 200°, the opposite is observed, with larger v_R for higher density regions of this azimuthal angle (R_phot ≲ 4.5 kpc).

Fig. 19 Results of the Fourier transform of stellar density maps of the LMC for the BL and the combined samples (upper and bottom panels, respectively). Left and middle panels: amplitude and phases of the first- and second-order asymmetries (dashed black and solid green lines, respectively). The axisymmetric density is shown as a dotted blue line. Right panel: k = 1 perturbation over the bisymmetry strength ratio. The galactocentric radius is given here in the photometric frame, with Ωphot = 10° as the reference photometric position angle.

Fig. 20 Stellar density map of the combined sample of the LMC. The lopsided eastern spiral arm is located beyond Rphot = 3.8 kpc (shownas an ellipse). The perpendicular axes show the orientation of the photometric major and minor axes (Sect. 7.2). The density scale shown by the colour bar is logarithmic; the values indicate the star count in pixels.

8 Conclusions

Using the new Gaia EDR3 data, we studied the structure and kinematics of the Magellanic Clouds with a new basis. The increased completeness and precision of the new release have allowed us to improve upon previous results using Gaia DR2, although (by design, because this is just a demonstration paper) we have certainly barely scratched the potential of the new data for the study of the Clouds.

In Sect. 3 we compared the Gaia DR2 and Gaia EDR3 data in the region of the LMC and SMC, showing the improvement in the astrometry and photometry from one release to the other. Not only the precision has increased, but the systematic effects are significantly reduced. The reduced crowding effects in the photometry are particularly relevant for the central regions of the Clouds.

We have explored the use of the astrometric data for the determination of the 3D structure of the LMC. Our attempts to use Bayesian modelling to reconstruct the geometry of the system have shown that despite the significant improvements in Gaia EDR3, the systematic effects still present on the parallaxes (the regional zero-points discussed in Lindegren et al. 2021b) distort mostof the signal of the 3D structure of the LMC (in the astrometry), and that there is not enough information in the summary statistics used by our approximate Bayesian method to simultaneously infer the local parallax zero-point variations and the geometry of the LMC. However, we do not rule out that it may be possible to determine this geometry by adding additional external restrictions and/or finding an optimal way to include additional information of the density distribution of the stars in the LMC area.

Our kinematic modelling of the proper motions has allowed us to derive radial and tangential velocity maps and global profiles for the LMC. This is the first time to our knowledge that the two planar components of the ordered and random motions are derived for multiple stellar evolutionary phases in a galactic disc outside the Milky Way. We show that younger stellar phases rotate fasterthan older ones. This is a clear effect of the asymmetric drift on the stellar kinematics. We have also been able to find the rotational centre of the stars in the LMC and showed that is significantly offset from the rotational centre of the H I gas.

On the other hand, we showed that the radial velocity is mostly negative in the inner 5 kpc, and inward radial motions are stronger in the bar region for younger stars. This velocity varies strongly as a function position with respect to the bar, and therefore to some extent reflects streaming motions along the bar. We observed asymmetric radial and tangential motions in the disc, such as those across and along the large-scale outer spiral arm of the LMC, in which the tangential velocity is larger in higher stellar density regions of the arm, and lower at lower density. There appears to be no clear rule in the streaming of the radial motion as a function of the position in the arm, however.

We showed that the radial and tangential random motions decrease from the disc centre out to the outskirts, but not at the same rate in the evolutionary phases. Older stars lie in a dynamically hotter disc than younger stars. While we have found higher velocity dispersions aligned with the bar in all stellar phases for both components, we also found evidence that only the radial component of more evolved stars exhibits a central feature of higher amplitude, which is oriented perpendicular to the bar.

Our analysis of the stellar density maps has shown more concentrated and clumpier distributions for younger stars in the bar and inner spiral structure(s) than older disc stars. Analysis of Fourier harmonics of density maps also revealed that the inner disc is perturbed predominantly by bisymmetries, the bar and spiral arms, while the outer disc is perturbed by a dominant lopsided spiral arm. The peak strengths of these perturbations can be as high as 40–60% of the axisymmetric amplitude in the inner region, and higher in the disc outskirts for the lopsidedness.

Finally, we carried out a specific study of the Magellanic Bridge. Using two different evolutionary phases (the young and RC samples), we were able to trace the density and velocity flow of the stars from the SMC towards the LMC following the Bridge. We showed that it apparently wraps around the LMC, connecting with the young southern arm-like structure of the galaxy. The quality of the Gaia EDR3 proper motions also allowed us to confirm the bridge in the RC evolutionary phase, at a position slightly shifted from that of the young evolutionary phase. Additionally, we were able to study the outskirts of both Magellanic Clouds and detected some well-known features, such as the north and south tidal arms of the LMC and the northern overdensity of the SMC. Our data also suggests a faint overdensity east of the LMC, which has only recently been reported with the help of near-infrared maps (El Youssoufi, priv. comm.).

Fig. 21 Streaming motions across and along the eastern spiral arm of the LMC. The azimuth-radius diagrams show the stellar density (top panel), the tangential velocity (middle panel), and radial velocity (bottom panel). Contours represent the stellar density shown in the upper panel.

Acknowledgements

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. This work presents results from the European Space Agency (ESA) space mission Gaia. Gaia data are being processed by the Gaia Data Processing and Analysis Consortium (DPAC). Funding for the DPAC is provided by national institutions, in particular the institutions participating in the Gaia MultiLateral Agreement (MLA). The Gaia mission website is https://www.cosmos.esa.int/gaia. The Gaia archive website is https://archives.esac.esa.int/gaia. The Gaia mission and data processing have financially been supported by, in alphabetical order by country: the Algerian Centre de Recherche en Astronomie, Astrophysique et Géophysique of Bouzareah Observatory; the Austrian Fonds zur Förderung der wissenschaftlichen Forschung (FWF) Hertha Firnberg Programme through grants T359, P20046, and P23737; the BELgian federal Science Policy Office (BELSPO) through various PROgramme de Développement d’Expériences scientifiques (PRODEX) grants and the Polish Academy of Sciences - Fonds Wetenschappelijk Onderzoek through grant VS.091.16N, and the Fonds de la Recherche Scientifique (FNRS); the Brazil-France exchange programmes Fundaçãode Amparo à Pesquisa do Estado de São Paulo (FAPESP) and Coordenação de Aperfeicoamento de Pessoal de Nível Superior (CAPES) - Comité Français d’Evaluation de la Coopération Universitaire et Scientifique avec le Brésil (COFECUB); the National Science Foundation of China (NSFC) through grants 11573054 and 11703065 and the China Scholarship Council through grant 201806040200; the Tenure Track Pilot Programme of the Croatian Science Foundation and the École Polytechnique Fédérale de Lausanne and the project TTP-2018-07-1171 “Mining the Variable Sky”, with the funds of the Croatian-Swiss Research Programme; the Czech-Republic Ministry of Education, Youth, and Sports through grant LG 15010 and INTER-EXCELLENCE grant LTAUSA18093, and the Czech Space Office through ESA PECS contract 98058; the Danish Ministry of Science; the Estonian Ministry of Education and Research through grant IUT40-1; the European Commission’s Sixth Framework Programme through the European Leadership in Space Astrometry (ELSA) Marie Curie Research Training Network (MRTN-CT-2006-033481), through Marie Curie project PIOF-GA-2009-255267 (Space AsteroSeismology & RR Lyrae stars, SAS-RRL), and through a Marie Curie Transfer-of-Knowledge (ToK) fellowship (MTKD-CT-2004-014188); the European Commission’s Seventh Framework Programme through grant FP7-606740 (FP7-SPACE-2013-1) for the Gaia European Network for Improved data User Services (GENIUS) and through grant 264895 for the Gaia Research for European Astronomy Training (GREAT-ITN) network; the European Research Council (ERC) through grants 320360 and 647208 and through the European Union’s Horizon 2020 research and innovation and excellent science programmes through Marie Skłodowska-Curie grant 745617 as well as grants 670519 (Mixing and Angular Momentum tranSport of massIvE stars – MAMSIE), 687378 (Small Bodies: Near and Far), 682115 (Using the Magellanic Clouds to Understand the Interaction of Galaxies), and 695099 (A sub-percent distance scale from binaries and Cepheids – CepBin); the European Science Foundation (ESF), in the framework of the Gaia Research for European Astronomy Training Research Network Programme (GREAT-ESF); the European Space Agency (ESA) in the framework of the Gaia project, through the Plan for European Cooperating States (PECS) programme through grants for Slovenia, through contracts C98090 and 4000106398/12/NL/KML for Hungary, and through contract 4000115263/15/NL/IB for Germany; the Academy of Finland and the Magnus Ehrnrooth Foundation; the French Centre National d’Etudes Spatiales (CNES), the Agence Nationale de la Recherche (ANR) through grant ANR-10-IDEX-0001-02 for the “Investissements d’avenir” programme, through grant ANR-15-CE31-0007 for project “Modelling the Milky Way in the Gaia era” (MOD4Gaia), through grant ANR-14-CE33-0014-01 for project “The Milky Way disc formation in the Gaia era” (ARCHEOGAL), and through grant ANR-15-CE31-0012-01 for project “Unlocking the potential of Cepheids as primary distance calibrators” (UnlockCepheids), the Centre National de la Recherche Scientifique (CNRS) and its SNO Gaia of the Institut des Sciences de l’Univers (INSU), the “Action Fédératrice Gaia”’ of the Observatoire de Paris, the Région de Franche-Comté, and the Programme National de Gravitation, Références, Astronomie, et Métrologie (GRAM) of CNRS/INSU with the Institut National Polytechnique (INP) and the Institut National de Physique nucléaire et de Physique des Particules (IN2P3) co-funded by CNES; the German Aerospace Agency (Deutsches Zentrum für Luft- und Raumfahrt e.V., DLR) through grants 50QG0501, 50QG0601, 50QG0602, 50QG0701, 50QG0901, 50QG1001, 50QG1101, 50QG1401, 50QG1402, 50QG1403, 50QG1404, and 50QG1904 and the Centre for Information Services and High Performance Computing (ZIH) at the Technische Universität (TU) Dresden for generous allocations of computer time; the Hungarian Academy of Sciences through the Lendület Programme grants LP2014-17 and LP2018-7 and through the Premium Postdoctoral Research Programme (L. Molnár), and the Hungarian National Research, Development, and Innovation Office (NKFIH) through grant KH_18-130405; the Science Foundation Ireland (SFI) through a Royal Society - SFI University Research Fellowship (M. Fraser); the Israel Science Foundation (ISF) through grant 848/16; the Agenzia Spaziale Italiana (ASI) through contracts I/037/08/0, I/058/10/0, 2014-025-R.0, 2014-025-R.1.2015, and 2018-24-HH.0 to the Italian Istituto Nazionale di Astrofisica (INAF), contract 2014-049-R.0/1/2 to INAF for the Space Science Data Centre (SSDC, formerly known as the ASI Science Data Center, ASDC), contracts I/008/10/0, 2013/030/I.0, 2013-030-I.0.1-2015, and 2016-17-I.0 to the Aerospace Logistics Technology Engineering Company (ALTEC S.p.A.), INAF, and the Italian Ministry of Education, University, and Research (Ministero dell’Istruzione, dell’Università e della Ricerca) through the Premiale project “MIning The Cosmos Big Data and Innovative Italian Technology for Frontier Astrophysics and Cosmology” (MITiC); the Netherlands Organisation for Scientific Research (NWO) through grant NWO-M-614.061.414, through a VICI grant (A. Helmi), and through a Spinoza prize (A. Helmi), and the Netherlands Research School for Astronomy (NOVA); the Polish National Science Centre through HARMONIA grant 2018/06/M/ST9/00311, DAINA grant 2017/27/L/ST9/03221, and PRELUDIUM grant 2017/25/N/ST9/01253, and the Ministry of Science and Higher Education (MNiSW) through grant DIR/WK/2018/12; the Portugese Fundação para a Ciência e a Tecnologia (FCT) through grants SFRH/BPD/74697/2010 and SFRH/BD/128840/2017 and the Strategic Programme UID/FIS/00099/2019 for CENTRA; the Slovenian Research Agency through grant P1-0188; the Spanish Ministry of Economy (MINECO/FEDER, UE) through grants ESP2016-80079-C2-1-R, ESP2016-80079-C2-2-R, RTI2018-095076-B-C21, RTI2018-095076-B-C22, BES-2016-078499, and BES-2017-083126 and the Juan de la Cierva formación 2015 grant FJCI-2015-2671, the Spanish Ministry of Education, Culture, and Sports through grant FPU16/03827, the Spanish Ministry of Science and Innovation (MICINN) through grant AYA2017-89841P for project “Estudio de las propiedades de los fósiles estelares en el entorno del Grupo Local” and through grant TIN2015-65316-P for project “Computación de Altas Prestaciones VII”, the Severo Ochoa Centre of Excellence Programme of the Spanish Government through grant SEV2015-0493, the Institute of Cosmos Sciences University of Barcelona (ICCUB, Unidad de Excelencia “María de Maeztu”) through grants MDM-2014-0369 and CEX2019-000918-M, the University of Barcelona’s official doctoral programme for the development of an R+D+i project through an Ajuts de Personal Investigador en Formació (APIF) grant, the Spanish Virtual Observatory through project AyA2017-84089, the Galician Regional Government, Xunta de Galicia, through grants ED431B-2018/42 and ED481A-2019/155, support received from the Centro de Investigación en Tecnologías de la Información y las Comunicaciones (CITIC) funded by the Xunta de Galicia, the Xunta de Galicia and the Centros Singulares de Investigación de Galicia for the period 2016-2019 through CITIC, the European Union through the European Regional Development Fund (ERDF) / Fondo Europeo de Desenvolvemento Rexional (FEDER) for the Galicia 2014-2020 Programme through grant ED431G-2019/01, the Red Española de Supercomputación (RES) computer resources at MareNostrum, the Barcelona Supercomputing Centre - Centro Nacional de Supercomputación (BSC-CNS) through activities AECT-2016-1-0006, AECT-2016-2-0013, AECT-2016-3-0011, and AECT-2017-1-0020, the Departament d’Innovació, Universitats i Empresa de la Generalitat de Catalunya through grant 2014-SGR-1051 for project ’Models de Programació i Entorns d’Execució Parallels’ (MPEXPAR), and Ramon y Cajal Fellowship RYC2018-025968-I; the Swedish National Space Agency (SNSA/Rymdstyrelsen); the Swiss State Secretariat for Education, Research, and Innovation through the ESA PRODEX programme, the Mesures d’Accompagnement, the Swiss Activités Nationales Complémentaires, and the Swiss National Science Foundation; the United Kingdom Particle Physics and Astronomy Research Council (PPARC), the United Kingdom Science and Technology Facilities Council (STFC), and the United Kingdom Space Agency (UKSA) through the following grants to the University of Bristol, the University of Cambridge, the University of Edinburgh, the University of Leicester, the Mullard Space Sciences Laboratory of University College London, and the United Kingdom Rutherford Appleton Laboratory (RAL): PP/D006511/1, PP/D006546/1, PP/D006570/1, ST/I000852/1, ST/J005045/1, ST/K00056X/1, ST/K000209/1, ST/K000756/1, ST/L006561/1, ST/N000595/1, ST/N000641/1, ST/N000978/1, ST/N001117/1, ST/S000089/1, ST/S000976/1, ST/S001123/1, ST/S001948/1, ST/S002103/1, and ST/V000969/1.

Appendix A Sky plots of the samples used

This appendix provides some additional figures complementary to the main text, collected here to avoid cluttering the main body and facilitate reading of the paper.

First, inFig. A.1 we illustrate the joint sky distribution of our two (LMC and SMC) basic clean samples. The selection radius for both clouds is clearly visible. Following this map, the next figures show the distribution of the mean G magnitude (Fig. A.2), mean G_BP − G_RP (Fig. A.3), the variation of this mean between Gaia DR2 and Gaia EDR3 (Fig. A.4), and the mean parallax (Fig. A.5).

Fig. A.1 Source density in the two circular areas with a radius of 20° and 11°, centred on the LMC and SMC. Only the 12.4M sources selected as potential members are included, and the criteria are therefore slightly different for the two circles.

Fig. A.2 Mean G magnitude in the LMC and SMC.

Fig. A.3 Mean colour, GBP − GRP in the LMC and SMC.

Fig. A.4 Mean change in colour, GBP − GRP, in Gaia DR2 and Gaia EDR3 for sources in the LMC and SMC. Positive values mean that the sources are now redder.

Fig. A.5 Mean parallax for sources in the LMC and SMC. No corrections have been applied to the parallaxes.

Next, we provide plots of the sky distribution for each one of our evolutionary phase samples; the plots show that each sample traces different structures of the clouds. Figure A.6 contains the plots for the LMC and Fig. A.7 those for the SMC.

Fig. A.6 Skydensity plots for the complete LMC clean sample and the different evolutionary phase subsamples.

Fig. A.7 Skydensity plots for the complete SMC clean sample and the different evolutionary phase subsamples.

Appendix B LMC velocity maps and profiles

We include here several figures that complement the analysis of the kinematics of the LMC presented in Sect. 5. Figures B.1–B.4 present the velocity maps of the LMC (azimuthal and radial velocities, mean values, and dispersion) for the eight evolutionary phases and the combined sample of stars. Finally, Fig. B.5 presents the velocity profiles of the LMC (as a function of the distance to its centre) traced using the different populations defined by our evolutionary phases.

Fig. B.1 Azimuthal velocity maps vϕ of the LMC for the combined sample (main panel) and the various evolutionary phases. The linear velocity scale shown by the colour bar in the main panel is the same in all subpanels and has been chosen to highlight velocity patterns better.

Fig. B.2 Radial velocity maps vR of the LMC for the combined sample (main panel) and the various evolutionary phases. The linear velocity scale shown by the colour bar in the main panel is the same in all subpanels and has been chosen to highlight velocity patterns better.

Fig. B.3 Azimuthal velocity dispersion maps σϕ of the LMC for the combined sample (main panel) and the various evolutionary phases. The linear velocity scale shown by the colour bar in the main panel is the same in all subpanels and has been chosen to highlight velocity patterns better.

Fig. B.4 Radial velocity dispersion maps σR of the LMC for the combined sample (main panel) and the various evolutionary phases. The linear velocity scale shown by the colour bar in the main panel is the same in all subpanels and has been chosen to highlight velocity patterns better.

Fig. B.5 Velocity profiles of the LMC. Rotation curve (vϕ), tangential velocity dispersion (σϕ), radial velocity (vR), and radial velocity dispersion (σR) from leftto right for each stellar evolutionary phase (from top to bottom). Velocities are in km s−1. Bottom row: result for the combined sample. The blue shaded areas correspond to the uncertainties.

References

All Tables

All Figures

	Fig. 1 Sky density plots for the LMC (left) and SMC (right) clean samples (after spatial and proper motion selection). Top row: plots in equatorial coordinates. Bottom row: orthographic projection (as used in Sect. 3)
In the text

	Fig. 2 Areas (as defined by the polygons given in the text) of the CMD for the LMC (left) and SMC (right) evolutionary phases. The colours are not corrected for reddening for the selection.
In the text

	Fig. 3 Estimated age distribution of the selected evolutionary phase. Based on a simulation using a constant star formation rate, the age-metallicity relation by Harris & Zaritsky (2009), and PARSEC1.2 models
In the text

	Fig. 4 Comparison of the parallaxes (left) and proper motions in the x and y directions (middle and right, respectively) of LMC sources in Gaia DR2 (upper panels) and Gaia EDR3 (lower panels).
In the text

	Fig. 5 Sameas in Fig. 4, but for the SMC.
In the text

	Fig. 6 Comparison of the residual proper motion fields of the LMC after a first-order approximation of the field was subtracted to emphasise the systematic errors in Gaia DR2 (upper panels) and Gaia EDR3 (lower panels).
In the text

	Fig. 7 Sameas in Fig. 6, but for the SMC.
In the text

	Fig. 8 Illustration of the proper motion field of the LMC (left) and SMC (right) using line-integral convolution. We set the alpha parameter (opacity) of the coloured lines according to the density, with the densest regions being the most opaque.
In the text

	Fig. 9 Photometric excess factor, i.e. the sum of fluxes in the GBP and GRP bands over the G flux. Left: for Gaia DR2. Right: for Gaia EDR3.
In the text

	Fig. 10 Positions of the centre of the LMC as found by different studies (as described in the text, SkyMapper estimates are referred to as “Sky” in the figure) superimposed on an image coloured according to the total density in star counts in the inner few degrees of the LMC. The centres found in this study lie closer to the photometric centre than to the centre of H I rotation.
In the text

	Fig. 11 Plots of vϕ and vR for our main model (top) and for our model with the centre left as a free parameter (bottom).
In the text

	Fig. 12 Example of velocity maps for the LMC. Left and right columns: young and evolved evolutionary phases (Young 2 and RC, respectively, see Sect. 2.3). Upper and bottom rows: vϕ and vR, respectively.The linear velocity scales shown by colour bars are the same for the two stellar evolutionary phases and were chosen to show the structure inside the velocity fields more clearly.
In the text

	Fig. 13 Same as in Fig. 12, but for the velocity dispersion.
In the text

Fig. 14 Stellar velocity curves of the LMC. Left and right panels: rotation curves and radial motions, respectively. Coloured lines are for the eight evolutionary phases and the combined sample, as given by the legend. In the left panel, blue circles, red diamonds and magenta squares are the rotation velocities for the samples of youngerand older stars from line-of-sight spectroscopic and HST astrometric measurements published in van der Marel & Kallivayalil (2014). The dashed blue lines are the tangential and radial velocity profiles derived by Wan et al. (2020) from Gaia DR2 astrometry.

In the text

	Fig. 15 Sky density plot for the Bridge region using the Young1 and Young2 evolutionary phases (left panel) and the RC evolutionary phase (right panel). The bin size is 0.06 deg in x and y. The colour bar is in log scale and the black contours are at the levels 0.1, 1, 5, and 70 (young evolutionary phases) and 0.25, 1.5, 15, and 165 (RC sample).
In the text

Fig. 16 Vector field of the proper motions in the Magellanic Clouds using theYoung 1 +Young 2 (left panel) and RC (right panel) samples. The coordinates are centred on the SMC. In the background, to guide the eye of the reader, we show the density in logarithmic scale. The dashed line in the left panel shows the location of the young bridge from Belokurov et al. (2017). The velocity vector field is only shown for bins with more than 10 (200) stars in theYoung 1 +Young 2 (RC) sample.

In the text

	Fig. 17 Number count map for the LMCout sample (left): the names of the most notable substructures are given and, for a better visualisation, the inner parts of the LMC and SMC were omitted. Total velocity dispersion map (right): the velocity vector field is shown as orange arrows.
In the text

Fig. 18 Differential density maps for different stellar populations in the LMC. Left column: density map in the BL evolutionary phase. Middle panels from top to bottom: density maps of the young phase, i.e. stars with age <400 Myr, Young 3 phase, i.e. main- sequence stars with age 1−1.5 Gyr, and oldphase, i.e. stars with age >1 Gyr. Right column: normalised difference between the left and middle panel maps, shown with a logarithmic stretch.

In the text

Fig. 19 Results of the Fourier transform of stellar density maps of the LMC for the BL and the combined samples (upper and bottom panels, respectively). Left and middle panels: amplitude and phases of the first- and second-order asymmetries (dashed black and solid green lines, respectively). The axisymmetric density is shown as a dotted blue line. Right panel: k = 1 perturbation over the bisymmetry strength ratio. The galactocentric radius is given here in the photometric frame, with Ωphot = 10° as the reference photometric position angle.

In the text

	Fig. 20 Stellar density map of the combined sample of the LMC. The lopsided eastern spiral arm is located beyond Rphot = 3.8 kpc (shownas an ellipse). The perpendicular axes show the orientation of the photometric major and minor axes (Sect. 7.2). The density scale shown by the colour bar is logarithmic; the values indicate the star count in pixels.
In the text

	Fig. 21 Streaming motions across and along the eastern spiral arm of the LMC. The azimuth-radius diagrams show the stellar density (top panel), the tangential velocity (middle panel), and radial velocity (bottom panel). Contours represent the stellar density shown in the upper panel.
In the text

	Fig. A.1 Source density in the two circular areas with a radius of 20° and 11°, centred on the LMC and SMC. Only the 12.4M sources selected as potential members are included, and the criteria are therefore slightly different for the two circles.
In the text

	Fig. A.2 Mean G magnitude in the LMC and SMC.
In the text

	Fig. A.3 Mean colour, GBP − GRP in the LMC and SMC.
In the text

	Fig. A.4 Mean change in colour, GBP − GRP, in Gaia DR2 and Gaia EDR3 for sources in the LMC and SMC. Positive values mean that the sources are now redder.
In the text

	Fig. A.5 Mean parallax for sources in the LMC and SMC. No corrections have been applied to the parallaxes.
In the text

	Fig. A.6 Skydensity plots for the complete LMC clean sample and the different evolutionary phase subsamples.
In the text

	Fig. A.7 Skydensity plots for the complete SMC clean sample and the different evolutionary phase subsamples.
In the text

	Fig. B.1 Azimuthal velocity maps vϕ of the LMC for the combined sample (main panel) and the various evolutionary phases. The linear velocity scale shown by the colour bar in the main panel is the same in all subpanels and has been chosen to highlight velocity patterns better.
In the text

	Fig. B.2 Radial velocity maps vR of the LMC for the combined sample (main panel) and the various evolutionary phases. The linear velocity scale shown by the colour bar in the main panel is the same in all subpanels and has been chosen to highlight velocity patterns better.
In the text

	Fig. B.3 Azimuthal velocity dispersion maps σϕ of the LMC for the combined sample (main panel) and the various evolutionary phases. The linear velocity scale shown by the colour bar in the main panel is the same in all subpanels and has been chosen to highlight velocity patterns better.
In the text

	Fig. B.4 Radial velocity dispersion maps σR of the LMC for the combined sample (main panel) and the various evolutionary phases. The linear velocity scale shown by the colour bar in the main panel is the same in all subpanels and has been chosen to highlight velocity patterns better.
In the text

	Fig. B.5 Velocity profiles of the LMC. Rotation curve (vϕ), tangential velocity dispersion (σϕ), radial velocity (vR), and radial velocity dispersion (σR) from leftto right for each stellar evolutionary phase (from top to bottom). Velocities are in km s−1. Bottom row: result for the combined sample. The blue shaded areas correspond to the uncertainties.
In the text