Revealing the Structure and Internal Rotation of the Sagittarius Dwarf Spheroidal Galaxy with Gaia and Machine Learning

We assembled two different catalogs that will be used later to train prediction models to determine D and v_los for all stars from the gaia_source table. The v_los training table consists of 2427 stars with observed line-of-sight velocities, v_los,true, either from gaia_source, the APOGEE survey, or one of the following publicly available catalogs: Ibata et al. (1997), Giuffrida et al. (2010), Frinchaboy et al. (2012), and McDonald et al. (2012). Common stars between the catalogs show consistent measurements in most cases. For stars with less than 10 km s⁻¹ of standard deviation between catalogs, we adopted the error-weighted average of their v_{los, true} as the final value. Stars with larger standard deviation values were rejected.

We composed the distance training table from more than 4000 RR Lyrae stars around the Sgr center from the SOS and VariClassifier catalogs published by DPAC making use of Gaia DR2 (Holl et al. 2018; Clementini et al. 2019; Rimoldini et al. 2019). Contaminants were removed following the method described in Rimoldini et al. (2019). Their observed distances, D_true, are derived from their pulsation modes, therefore making them independent from the astrometric properties of the star. We estimated an average error of 1.16 kpc in D_true from M54 member stars that accounts for both observational and calibration effects (see Appendix C.1.1). Details on how the joint SOS and VariClassifier catalog was assembled and the computation of the RR Lyrae distances are available in Mateu et al. (2020).

Undesired systematic bias effects can occur when the training sample does not fully cover the entire space defined by x from which later predictions will be made. To avoid such problems, we initially considered a wider PM and ϖ range, as well as a wider region in the sky for the training lists than for the main catalog. Specifically, we selected all stars within 77 of the central coordinates of Sgr with PMs and ϖ compatible with Sgr by ±1 mas yr⁻¹ and ±0.5 plus 3σ mas, respectively. We then cleaned the two catalogs of contaminants and poorly measured stars following the exact same criteria as for the gaia_source data, ensuring that the three catalogs are subsamples representative of the same set of stars. Up to 3422 RR Lyrae and 2310 stars remained in the distance and line-of-sight velocity training sets, respectively, after the first membership selection. As a last step before fitting the models to the training lists, we standardized all variables in x by subtracting their averages and dividing them by their standard deviation.

The two training sets were used to create models that predict distances and line-of-sight velocities on the gaia_source table.

To obtain D, we used a k-neighbors regressor (KNR) model. The KNR models predict y_pred by local interpolation of the y_true of the k nearest neighbors in the parameter space defined by x in the training set. This interpolation is, in our case, weighted by the inverse of the euclidean distance to the target star in the parameter space defined by x . The KNR models are highly scalable and perform fast on large and well-distributed samples. The procedure is as follows.

The KNR model is fitted to the training distances, y_true, using sky coordinates and a set of astrometric magnitudes as x . The number of neighbors, k, and best combination of independent variables in x are determined through a grid search within a nested cross-validation (NCV) of the training set minimizing (a) the variance in the differences between the predicted and the true values for the distance (generalization error) and (b) possible artificial gradients along the sky introduced by the model (see Appendix C). The validated model is then used to predict distances on the gaia_source based on the same independent variables.

Choosing a KNR over other models has the advantage that the model regularizes (flattens) y_true when k > 1 is used, avoiding overfitting the data. The relatively large errors observed in our distances (1.16 kpc) and the large FWHM of their distribution (3.69 kpc) compared to previous works (Hamanowicz et al. 2016) indicate that regularization and averaging of the training sample are required prior to deriving distances for the main table (for more information about the errors in the training list and the KNR performance, see Appendix C.1).

The procedure to derive v_los is similar, but due to the sparse nature of the data (their distribution in the sky is far from uniform) and their nonlinearity, we decided to use a stacked regressor (SR; Wolpert 1992; Breiman 1996). The SR models take advantage of different learners by stacking their output through a final regressor that computes the final prediction. They perform much slower than each of the estimators separately, but their results are at least as good as the best of the stacked estimators while improving the results in most cases. We combined six learners, three nonlinear (extra trees, support vector machine with a radial basis function, and Gaussian process) and three linear (lasso, support vector machine with a linear kernel, and elastic net), and used a neural network to combine their outputs. The SR was trained through an out-of-fold scheme: in k = 5 successive iterations, k-1 folds are used to fit the first layer of regressors to make predictions on the remaining subset. The predictions are then stacked and provided as input to the final metaregressor. The optimal metaparameters of the six learners inside our SR model, as well as the best combination of independent variables in x , are optimized through a grid search over the metaparameter space within an NCV. The resulting model was able to reproduce nonlinear features while showing stability against high-order oscillations in areas where no information is available, such as the outer regions of the galaxy, where the data density is much lower. (For more information about the SR architecture and performance, see Appendix C.2.)

The procedure to derive y_pred and its associated error is done in a Monte Carlo fashion, considering and propagating the errors from all independent and dependent variables. A total of 10⁶ realizations were performed, randomly sampling each variable from a normal distribution centered on its corresponding nominal value. The predicted values using the nominal values of the independent variables are assumed to be the nominal values of the final result. The quadrature addition of the standard deviation of all realizations plus the square root of the generalization error derived from the NCV test is adopted as the total error.

Stars in the main catalog and the two training sets are finally selected based on their full velocity space, i.e., 3D velocities, through a 3D Gaussian using the same rejection criteria adopted in Section 2. No rejection was imposed on the spatial coordinates. The whole procedure is repeated until no more stars are rejected and the solution converges. The final lists contain 3305, 1992, and 120,168 stars for the RR Lyrae, the v_los, and the final clean sample, respectively.

Values for k and the best combination of independent variables, x , vary a bit from one iteration to another for the KNR. A larger k produces more stable results but at the cost of flattening the results too much, artificially removing gradients from the data. Our criteria seem to be a good compromise, minimizing such changes in distance gradients and keeping low deviations from the true values. The best results were obtained using k between 10 and 30, and x = (x, y, μ_x , μ_y , ϖ), where x, y, μ_x , and μ_y are the orthographic projections of the usual celestial coordinates and PMs (see Appendix C.1.2). The optimal metaparameters of the SR model, as well as the best x , also vary between iterations yet produce fully compatible results. See Appendix C for more information.

After convergence, we evaluated the performance of both models by comparing y_true and y_pred from the NCV tests on each corresponding final training list. The NCV procedure uses a series of train, validation, and test set splits to effectively optimize and test the model over different data sets, thus allowing an independent assessment of the models' ability to predict on new data.

Figure 2 shows the predicted distances, D, for our final clean sample and the residuals from the NCV test, D_true – D_pred, on the coordinate system (Λ, B), defined along the stream (Law & Majewski 2010). Here Λ almost coincides with the optical major axis of the system (as projected on the sky) and B with the minor one (slightly offset with respect to the center of the galaxy). A linear fit to the training set, y_true, has been added to shown the general tendency: the western regions of the galaxy (negative Λ) are, on average, closer to the Sun than the eastern parts (positive Λ). A less obvious trend can also be seen along B, with negative B closer to us.

Zoom In Zoom Out Reset image size

The agreement between the predicted and true distances is very good, with both distributions showing the same average behavior and a constant variance between D_true and D_pred along both axes. Some departures from the linear model (red dashed line) can be easily spotted in our predicted distances for the clean sample (contours), especially along Λ. These departures are also present on the training sample and point to a real twist in the distance along the line of sight. The KNR does a good job regularizing the distances, with a line-of-sight thickness FWHM ∼ 2.2 kpc for horizontal-branch (HB) stars around the instability strip. This is closer to the 2.42 kpc found by the OGLE team (Hamanowicz et al. 2016) and our real FWHM estimation of 2.48 kpc than the FWHM measured directly over our RR Lyrae sample (3.69 kpc). These are not worrisome differences for our purposes of analyzing the general distance trends in Sgr, further taking into account that a number of factors could be affecting the measured thickness along the line of sight, such as the number of stars used, metallicity calibration, etc. See Appendix C.1.1 for more information.

Figure 3 shows the predicted and true v_los for the clean final sample and the training sample, respectively. Again, some departure from the linear fit can be observed in both samples, especially in the very central regions (∣Λ∣ < 1°) and the most external ones (∣Λ∣ > 4°). These changes roughly coincide with the points at which changes in distance are also present, indicating the presence of tidal tails with a slightly different kinematics than the core of the galaxy.

Zoom In Zoom Out Reset image size

Interestingly, the standard deviation of the differences between the predicted v_los and the real one was found to be σ(v_los) = 11.7 km s⁻¹, just above the internal velocity dispersion of the galaxy. This, combined with the small observational errors in the training sample and the fact that the spatial distribution of this quantity along the sky is very similar to the one from intrinsic dispersion, suggests that the model fails to reproduce the internal v_los dispersion of Sgr but reproduces well the general trends in the galaxy. This comes as no surprise, since the model has been regularized to avoid overfitting, which also limits its capacity to fit random velocities. We take that into account during our subsequent Monte Carlo experiments by allowing the stars to have a v_los in a range that includes this dispersion. For distances, differences were of the order of the error in distance for individual RR Lyraes. In this case, the largest differences were observed for the most extreme cases, i.e., stars with the largest or smallest distances in the training sample, indicating that it is a systematic effect caused by the regularization of the predicted distances, rather than random errors. More information can be found in Appendix C.2.

Despite all of our efforts to clean up our sample from nonmember stars, some MW contaminants may have passed through our screening process. To assess this contamination rate, we compared our member star list with an equivalent sample from the Besancon 2016 model of the Galaxy using the Gaia Object Generator (GOG) model (Luri et al. 2014). We simulated DR2 errors in the GOG model and performed the exact same selections as with real data (see Appendix A for further information). This resulted in ∼7500 GOG stars passing through our membership selection method, which indicates a ∼ 6% of possible contamination from MW stars in our final sample. Figure 4 shows the final distribution in the sky of the stars selected as members and rejected for the gaia_source table and GOG model.

Zoom In Zoom Out Reset image size

The stars from the GOG model that have passed our selection criteria are stars whose phase space, color, and magnitude are indistinguishable from Sgr member stars, not only in bulk but by sectors along the sky. Their contribution to our final measurements is well below the 1σ level of uncertainty of these; thus, we do not expect any critical impact from their possible presence in our results. It is worth noting that some probable Sgr members are rejected (overdensity in the center area of the bottom left panel in Figure 4). These are mostly stars with large astrometric and photometric errors. Their distribution in the sky follows the stellar density profile of Sgr, and their PMs randomly scatter in the vicinity of the cluster of members in Figure 6. These stars are expected to have a null contribution to the bulk dynamical properties of the galaxy given their small number and the fact that their astrometric properties are indistinguishable from those of Sgr in bulk and their surrounding neighbor stars within 1° (see Section 2).

The final selection of member stars is shown in the PM–parallax space in Figure 6. The set of stars has its center at (μ_α⋆, μ_δ , ϖ) = ( −2.6924 ± 0.0006, −1.3713 ± 0.0005, 0.0176 ± 0.0003), with some correlation between parallax and PMs, which makes the simultaneous fitting of the three variables necessary in order to maximize the member/nonmember ratio. Our algorithm does a good job of rejecting MW stars (∼2.7 × 10⁶) while keeping a high fraction of Sgr stars in the sample (∼1.4 × 10⁵).

Zoom In Zoom Out Reset image size

Our membership selection method iteratively finds the centroid of all stars within 2σ of a 6D Gaussian fitted to their galactocentric Cartesian phase space (X, Y, Z, v_X , v_Y , v_Z ) as the error-weighted average of such coordinates. The resulting centroid is assumed to be the COM of Sgr, and its coordinates are then transformed to sky coordinates (R.A., decl., D, μ_α⋆, μ_δ , v_los). Such a transformation is done assuming the distance of the Sun from the Galactic center and the circular velocity of the local standard of rest (LSR) to be R₀ = 8.29 ± 0.16 kpc and V₀ = 239 ± 5 km s⁻¹, respectively (McMillan 2011). The solar peculiar velocity with respect to the LSR was taken from the estimates of Schönrich et al. (2010): (U_pec, V_pec, W_pec) = (11.10, 12.24, 7.25) km s⁻¹ with uncertainties of (1.23, 2.05, 0.62) km s⁻¹.

The M54 stars could be shifting the centroid of Sgr toward the center of the globular cluster. In order to study this influence, we recalculated the COM for our table upon removal of all stars around M54 up to its tidal radius. We assumed r_t = r_c 10^c , where c = 2.04 and r_c = 009 are the concentration and core radius of the cluster (Harris 1996, 2010 edition). This provided a list of 1439 stars, from which we selected 352 very likely members of M54 based on their position on the CMD.

In Table 1, we list the sky coordinates and velocities for the COM of all member stars from the last iteration of our selection method. The first row shows the COM for the whole sample, including M54 stars. The second row shows the COM after removing M54 stars within its tidal radius. The third row shows the COM for M54 based on our 352 stars selected. The uncertainties listed here were obtained from a Monte Carlo scheme that propagates all observational uncertainties and their correlations, including those for the Sun and possible systematic errors (see Appendix D).

The presence of M54 stars seems to have a negligible influence on the derived COM properties of Sgr, except for its coordinates in the sky, (α, δ). Such a difference is not surprising; M54 stars are highly concentrated on the sky compared to the rest of the Sgr stars and therefore have a high weight when deriving the Sgr centroid in the sky.

The large difference in coordinates between Sgr and M54 (${\rm{\Delta }}\alpha =16\buildrel{\,\prime}\over{.} 9\pm 2\buildrel{\,\prime}\over{.} 2,{\rm{\Delta }}\delta =10\buildrel{\,\prime}\over{.} 8\pm 0\buildrel{\,\prime}\over{.} 6$) is also expected given the different areas considered. The systematics known to be affecting our sample produce a nonuniform surface stellar density along the observed region, thus affecting the centroid. If we focus on the rest of the properties, Sgr and M54 seem to be indistinguishable within the errors. Furthermore, all M54 properties are compatible with those derived for Sgr after removing M54 stars, indicating that both stellar populations share position and dynamics. We will, therefore, consider M54 as part of Sgr regarding astrometry in any respect.

We also studied the effect of large-scale systematics in the PMs of Sgr using quasars (see Appendix D.2). A total of 217 quasars were used to derive the median zero-points in PMs in the field of Sgr. We obtained $\tilde{{\mu }_{\alpha \star }}=0.03\pm 0.05$ and $\tilde{{\mu }_{\delta }}\,=0.09\pm 0.05\ \mathrm{mas}\ {\mathrm{yr}}^{-1}$, which indicate that Sgr's motion may be slower toward the north direction. The corrected PMs of Sgr and M54 are listed in Table 1 in parentheses.

Table 2 lists the galactocentric positions and velocities for Sgr. The core of Sgr is at a distance of 18.11 ± 0.15 kpc from the Galactic center, almost at the opposite side from the Sun location and 6.383 ± 0.007 kpc below the Galactic plane.

Currently, Sgr is on its way to cross the Galactic plane with a net velocity of 311.1 ± 1.1 km s⁻¹ (313 ± 4 if we correct from large-scale systematics) with respect to the Galactic center. Most of its velocity is tangential, V_T = 276.0 ± 1.2 km s⁻¹ (280 ± 6 km s⁻¹), while Sgr currently moves away from the Galactic center with a radial velocity of 143.6 ± 1.7 km s⁻¹ (140 ± 2). In general, its motion is well aligned with the Sgr stream, which lies close to the X–Z plane of the MW (misaligned by ∼ 15°). The velocity and position of Sgr indicate that it must have passed through its orbit pericenter very recently, and it is now receding with respect to the MW center.

The position of any point in Sgr can be uniquely determined by its R.A. and decl. on the sky, (α, δ), and its distance, D. To simplify our analysis and have a better interpretation of our results in terms of internal kinematics, we introduced a new reference system centered on the galaxy COM with coordinates (α₀, δ₀, D₀). Angular coordinates (ρ, ϕ) are defined on the celestial sphere, where ρ is the angular distance between the points (α, δ) and (α₀, δ₀) and ϕ is the position angle (PA) of the point (α, δ) with respect to (α₀, δ₀). The angle ϕ is measured counterclockwise starting from the axis that runs in the direction of decreasing R.A. at constant decl. δ₀ (see Figure 1 from van der Marel et al. 2002). Equations (1)–(3) from van der Marel & Cioni (2001) allow (ρ, ϕ) to be calculated from any (α, δ):

$$\begin{eqnarray}\begin{array}{c}\begin{array}{rcl}\rho & = & \arccos (\cos \delta \cos {\delta }_{0}\cos (\alpha -{\alpha }_{0})+\sin \delta \sin {\delta }_{0})\\ \phi & = & \arctan 2\left(\displaystyle \frac{\sin \delta \cos {\delta }_{0}-\cos \delta \sin {\delta }_{0}\cos (\alpha -{\alpha }_{0})}{-\cos \delta \sin (\alpha -{\alpha }_{0})}\right).\end{array}\end{array}\end{eqnarray} \tag{ 2 }$$

The velocity vector at any position (ρ, ϕ, D) can be decomposed into three orthogonal components:

$$\begin{eqnarray}&&{v}_{1}\equiv \displaystyle \frac{{dD}}{{dt}},{v}_{2}\equiv D\displaystyle \frac{d\rho }{{dt}},{v}_{3}\equiv D\sin \rho \displaystyle \frac{d\phi }{{dt}},\end{eqnarray} \tag{ 3 }$$

where v₁ is v_los and v₂ and v₃ are the radial and tangential components of the velocity with respect to the COM in the plane of the sky. These components can be derived in kilometers per second from any (μ_α*, μ_δ ) as

$$\begin{eqnarray}\begin{array}{rcl}{v}_{1} & = & {v}_{\mathrm{los}}\\ {v}_{2} & = & 4.74D[{\mu }_{\alpha * }\sin {\rm{\Gamma }}+{\mu }_{\delta }\cos {\rm{\Gamma }}]\\ {v}_{3} & = & 4.74D[{\mu }_{\alpha * }\cos {\rm{\Gamma }}-{\mu }_{\delta }\sin {\rm{\Gamma }}],\end{array}\end{eqnarray} \tag{ 4 }$$

where Γ is the rotation angle of (v₂, v₃) to the frame on the sky, which is determined by

$$\begin{eqnarray}\begin{array}{rcl}\sin {\rm{\Gamma }} & = & [\cos {\delta }_{0}\sin (\alpha -{\alpha }_{0})]/\sin \rho \\ \cos {\rm{\Gamma }} & = & [\sin \delta \cos {\delta }_{0}\cos (\alpha -{\alpha }_{0})-\cos \delta \sin {\delta }_{0}]/\sin \rho .\end{array}\end{eqnarray} \tag{ 5 }$$

The COM motion must be subtracted prior to the analysis of internal kinematics. Its contribution to the velocity field can be derived as

$$\begin{eqnarray}\begin{array}{rcl}{v}_{1} & = & {v}_{t}\sin \rho \cos (\phi -{\theta }_{t})+{v}_{\mathrm{los},0}\cos \rho \\ {v}_{2} & = & {v}_{t}\cos \rho \cos (\phi -{\theta }_{t})-{v}_{\mathrm{los},0}\sin \rho \\ {v}_{3} & = & -{v}_{t}\sin (\phi -{\theta }_{t}),\end{array}\end{eqnarray} \tag{ 6 }$$

where v_los,0 is the systemic v_los of the galaxy, v_t is the velocity of the COM projected on the sky, and θ_t is the angle that defines the direction of v_t defined using the same criterion as ϕ:

$$\begin{eqnarray}\begin{array}{rcl}{v}_{t} & = & 4.74{D}_{0}\sqrt{{\mu }_{\alpha * 0}^{2}+{\mu }_{\delta 0}^{2}}\\ {\theta }_{t} & = & \arctan \left(\displaystyle \frac{{\mu }_{\delta 0}}{-{\mu }_{\alpha * 0}}\right).\end{array}\end{eqnarray} \tag{ 7 }$$

We adopted all values listed in Table 1 in order to derive the COM contribution to v₁, v₂, and v₃.

The Cartesian coordinates (x, y, z) can be derived from (ρ, ϕ, D) as

$$\begin{eqnarray}\begin{array}{rcl}x & = & \sin \rho \cos \phi D\\ y & = & \sin \rho \sin \phi D,\\ z & = & {D}_{0}-\cos \rho D,\end{array}\end{eqnarray} \tag{ 8 }$$

while their respective differentials (v_x , v_y , v_z ) can be obtained from (v₁, v₂, v₃) as

$$\begin{eqnarray}\begin{array}{rcl}{v}_{x} & = & {v}_{1}\sin \rho \cos \phi +{v}_{2}\cos \rho \cos \phi -{v}_{3}\sin \phi \\ {v}_{y} & = & {v}_{1}\sin \rho \sin \phi +{v}_{2}\cos \rho \sin \phi +{v}_{3}\cos \phi .\\ {v}_{z} & = & -{v}_{1}\cos \rho +{v}_{2}\sin \rho .\end{array}\end{eqnarray} \tag{ 9 }$$

We derived (v₁, v₂, v₃) for all of the stars in our sample. As expected after the subtraction of the COM motion, the median of the line-of-sight component, $\tilde{{v}_{1}}=0.1\pm 0.1\ \mathrm{km}\ {{\rm{s}}}^{-1}$, ⁷ is compatible with zero. Figure 7 shows the distribution of v₂ and v₃ for all of the stars in the sample. The median radial component of the stellar velocities in the sky, $\tilde{{v}_{2}}\,=-1.1\pm 0.2\ \mathrm{km}\ {{\rm{s}}}^{-1}$, indicates that the galaxy is contracting, while its tangential component, $\tilde{{v}_{3}}=2.9\pm 0.2\ \mathrm{km}\ {{\rm{s}}}^{-1}$, indicates that Sgr is rotating counterclockwise as projected on the celestial sphere. The medians of the Monte Carlo realizations of the same quantities were closer to zero but still significant in rotation: ⁸ 0.1 ± 0.1, −0.15 ± 0.17, and 1.71 ± 0.16 km s⁻¹ for $\tilde{{v}_{1}}$, $\tilde{{v}_{2}}$, and $\tilde{{v}_{3}}$, respectively. The distribution of v₂ and v₃ extends toward positive and negative values of v₂ and v₃, respectively. This is likely to be produced by the influence of tidal tails in the kinematics of the galaxy. Indeed, stars with higher values of v₂ and lower values of v₃ are located at larger distances along the stream. Stars with v₂ > 0 km s⁻¹ and v₂ < −50 km s⁻¹ were located at a median absolute $\tilde{{\rm{\Lambda }}}=2\buildrel{\circ}\over{.} 65$, whereas the rest are closer to the Sgr center: median absolute stream longitude $\tilde{{\rm{\Lambda }}}=1\buildrel{\circ}\over{.} 96$.

Zoom In Zoom Out Reset image size

Velocity maps provide an immediate qualitative view of the dynamics of Sgr. Figure 8 shows the stellar density, as well as the error-weighted average of each velocity component, (v_x , v_y , v_z ), for stars within Voronoi cells of approximately 1000 stars defined on their sky positions. Voronoi tessellation ensures a constant number of stars per unit of resolution ("pixel"), thus providing a constant signal-to-noise ratio (Cappellari & Copin 2003).

Zoom In Zoom Out Reset image size

The stellar density projected in the sky has an elliptical profile whose concentration peaks at the position of M54. Results also show a clear velocity gradient (>10 km s⁻¹) in v_z to first order comparable with the line-of-sight velocity, indicating rotation about the projected optical minor axis of Sgr. This gradient is almost parallel to the bulk tangential velocity of Sgr and the direction of the stream (along Λ).

Counterclockwise rotation in the plane of the sky is also clearly visible in v_x and v_y , where velocity gradients are inclined with respect to the axes due to the PA of Sgr in the sky (1024). Southern regions of the galaxy are moving west (right), while northern regions are moving east (left). Here v_y shows the same pattern, with eastern regions moving south and western ones moving north. Interesting features can be noticed in the v_x and v_y maps, such as the apparent contraction of the core, the two negative bands (v_x < 0) coming out toward positive values of Λ (east direction), or the region with negative v_y around (x, y) = (1.5, 1) kpc. These are the consequences of the rotation of the galaxy, the inclination of its main body with respect to the plane of the sky, and the presence of tidal tails moving outward from its core. It is also worth noticing an area at (x, y) ∼ (−2, 0) kpc with unexpected positive values of v_y and negative v_x . This is a region showing unusual small values of μ_δ compared to the average caused by Gaia systematic errors, resulting in clockwise rotation (see Appendix D.2). Stripes with the same inclination can be seen in Figures 1 and 4 and are related to the scanning pattern of the Gaia satellite. This band would probably show inverted velocities if it were not affected by errors. Systematics also affect the stellar density profiles, although they tend to average out when considering large areas, not changing the general view of the galaxy.

7.2.1. Sky Projections

The full phase space allows us to represent Sgr in all dimensions (6D), providing a better understanding of the structure and dynamics of the galaxy. Figure 9 shows the projected stellar density and velocities for all stars in our sample in the coordinate system defined by Equations (8) and (9). In this representation, the galaxy is located inside a cube whose sides are the planes x–y, x–z, and z–y. The x–y plane is aligned with the sky, with the x-axis pointing to the west and the y-axis to the north (same as in Figure 8). In planes x–z and z–y, the galaxy is projected along the line of sight, with positive values of z pointing to the Sun. The observer is therefore located at (0, 0, ∞).

Zoom In Zoom Out Reset image size

Agglomerative Voronoi tessellation was applied to the data, creating cells of approximately 1000 stars for which the median values of the star velocities were calculated. These median velocities are represented by black arrows and qualitatively show how different parts of the galaxy move with respect to the COM. We also have calculated the 3D internal angular momentum of Sgr, L = I ω, where I is the inertia tensor and ω is the angular speed, for stars in the center of Sgr within 1 kpc radius. It is shown by a blue arrow.

The main body of Sgr is inclined with respect to the plane of the sky, being closer to the Sun for positive values of x and y. The data also suggest the presence of a bar-like structure more than 2 kpc long, from the tips of which depart the tidal tails that form the stream. This bar-and-tails configuration confers on Sgr the "S" shape observed in the x–z plane, which also projects on the x–y plane. We should note that the tails depart the main body at distances of 1.5 kpc (∼3° in the sky), which is much closer to the center than the half-light radius found in the literature for Sgr (57; McConnachie 2012). Studies analyzing the properties of the core of the galaxy should be aware of the effects of such tails.

Rotation is most obvious on the x–z plane (counterclockwise), with some rotation signal also visible on the x–y plane (counterclockwise). The stripe showing unusually high μ_δ is clearly visible in the x–y plane as arrows pointing in the northeast direction (opposite from their neighbors). The least obvious rotation is seen in the z–y plane, with most of the stars located at positive values of y moving toward positive values of z (right) and stars with negative y moving toward negative values of z (left). The apparent expansion observed in this plane is a projection effect; stars at positive z are moving toward positive z, while stars at negative z are moving toward negative z (see x–z panel). The angular momentum of the galaxy, projected in each plane by a blue arrow, reflects these general rotation patterns.

Combined, all dimensions depict an inclined system that is closer to us on the west and north. Closer areas are moving in our direction, while the furthest are moving even further away. Of course, this is in the COM frame, as the whole galaxy is moving away from the Sun at v_los = 142.7 ± 1.7 km s⁻¹, much larger than the ∼20 km s⁻¹ gradient discussed here.

7.2.2. Galactocentric Projections

The most important external factor affecting Sgr morphology and kinematics is the MW potential. Tidal forces stretch the body of the galaxy, accelerating its stars depending on their relative distance to the MW center. Figure 10 shows Sgr projected on the galactocentric frame, (X, Y, Z), with the observer located at a infinite distance above the MW plane in (16.73, 2.44, ∞) kpc. The X–Y plane is almost aligned to the orbital plane of Sgr.

Zoom In Zoom Out Reset image size

The largest component of the bulk velocity for Sgr is found in the X–Z plane, which is also the plane where the largest velocity gradient is observed. This suggests that tidal forces are greatly affecting Sgr kinematics, shearing its body and stripping its stars from outer regions. On the other hand, some internal rotation is apparent in the innermost parts of the body, especially in the X–Z and Z–Y planes.

7.2.3. Principal Axis Projections

We can better understand Sgr's inner structure and kinematics by aligning the system with its principal axes of inertia, $(x^{\prime} ,y^{\prime} ,z^{\prime} )$. We calculated $(x^{\prime} ,y^{\prime} ,z^{\prime} )$ within a 3D sphere of a given radius r_3D, as well as the Euler angles necessary to align Sgr from its orientation in the sky to these axes, α, β, and γ, defined as z − y − x rotation.

The ratio between axes changes depending on the considered sphere radius, as does the orientation. We derived the direction of the principal axes of inertia of Sgr and the ratio of the FWHM of the stellar density along them for stars within spheres of radii 0.25 kpc ≤ r_3D ≤ 3 kpc (notice that no star is located farther than ∼2.8 kpc in our data). This is shown in Figure 11.

Zoom In Zoom Out Reset image size

The Sgr is a triaxial system at all radii considered in our sample. However, the ratio between its intermediate and shortest axes, given by $y^{\prime}$ and $z^{\prime}$, is never larger than 1.14, with an average of 1.08, indicating that Sgr is almost prolate within 3 kpc. The inclination of the longest axis, measured by β, remains almost constant up to distances of 0.75 kpc ∼ r_3D, as expected, since the bar dominates the stellar distribution at small radii. The angle α measures the inclination of the sky-projected major axis of the galaxy from the x-axis (east direction) in the counterclockwise direction, i.e., α = 270° − PA. Its variation within just 1 kpc, ∼60°, is due to the transition from the core of the galaxy to the tidal tails that elongate the shape of the galaxy along the stream direction. Here α stabilizes at around r_3D = 2.5 kpc at 1675 (PA = 1024). Lastly, γ measures the pitch angle, i.e., the inclination in the north–south direction after the β rotation has been applied. The largest variation in γ occurs at 0.5 kpc ∼ r_3D ≲ 0.75 kpc, radii for which the galaxy changes from having its north regions closer to us to the opposite. Given the evident change in morphology for radii larger than 1 kpc, we decided to adopt this value to calculate the nominal orientation of the Sgr principal axes. Using such a radius, the principal axes of Sgr, $(x^{\prime} ,y^{\prime} ,z^{\prime} )$, form net angles of (43°, 2°, 46°) with respect to the plane of the sky and (2°, 18°, 72°) with respect to the orbital plane of Sgr around the MW. We should emphasize, though, that these angles vary depending on the considered galactocentric radius, r_3D.

Figure 12 shows the projected stellar density and kinematics over the principal axes of inertia of Sgr. The Euler angles necessary to align Sgr from its orientation in the sky to its principal axes of inertia are α = 177° ± 3°, β = −40° ± 2°, and γ = 16° ± 2°, measured at r_3D = 1 kpc. These rotations align Sgr's longest principal axis $x^{\prime}$ with the x-axis, the intermediate $y^{\prime}$ with the y-axis, and the shortest axis $z^{\prime}$ with the z-axis, with (x, y, z) defined by Equation (8).

Zoom In Zoom Out Reset image size

In this projection, the galaxy is aligned with the bar, clearly visible in the $x^{\prime}$–$y^{\prime}$ and $x^{\prime}$–$z^{\prime}$ planes. The almost spherical view of Sgr from its longest axis, i.e., projected over its intermediate and shortest axes, shows the little length difference found between the intermediate and shortest axis ($z^{\prime}$–$y^{\prime}$ plane). Indeed, the ratio between the FWHM of the stellar distribution along the three principal axes inside the 1 kpc radius sphere is 1:0.67:0.60, not far from prolate.

To further investigate the nature and strength of the bar, we measured the bar mode A₂ from the positions of the stars projected onto the $x^{\prime}$–$y^{\prime}$ plane (along the shortest axis). In general, the amplitude of the mth Fourier mode of the discrete distribution of stars in the galaxy is calculated as ${A}_{m}=(1/N)\left|{{\rm{\Sigma }}}_{j\,=\,1}^{N}\exp ({im}{\phi }_{j})\right|$, where ϕ_j are the particle phases in cylindrical coordinates and N is the total number of stars (Sellwood & Athanassoula 1986; Debattista & Sellwood 2000). Values larger than 0.2 of the m = 2 term of this expansion indicate the presence of a strong bar. With A₂ = 0.36, our data suggest that the elongated shape of Sgr is indeed a strong bar. We should also notice that we obtained A₂ = 0.3 measured directly on the positions in the sky, also larger than 0.2.

The Sgr rotates mainly around its intermediate principal axis, $y^{\prime}$. This is clearly visible from the direction of the angular momentum, L, almost parallel to $y^{\prime}$, and the projected velocity field on the $x^{\prime}$–$z^{\prime}$ plane. Counterclockwise rotation can be observed in the $x^{\prime}$–$y^{\prime}$, causing the angular momentum, L, to point toward positive values of $z^{\prime}$. The galaxy also shows very little clockwise rotation in its inner regions about its major axis, $x^{\prime}$ ($z^{\prime}$–$y^{\prime}$ plane). The expansion of the outer regions of Sgr is most visible in the $x^{\prime}$–$y^{\prime}$ plane, where areas dominated by the arms are being pulled away from the main body of the galaxy.

We computed the absolute net velocity along each one of the principal axes of Sgr as

$$\begin{eqnarray*}{v}_{\mathrm{rad},i}=\left\{\begin{array}{lc}{v}_{\mathrm{rad},i} & \mathrm{if\; i}\geqslant 0\\ -{v}_{\mathrm{rad},i} & \mathrm{if\; i}\lt 0\end{array}\right.,\end{eqnarray*}$$

with i being each of the principal axes $x^{\prime}$, $y^{\prime}$, or $z^{\prime}$. These velocities as a function of the absolute value of the radial distance along their respective axis ∣i∣ are shown in Figure 13. The rotation velocity around the same axes, computed over the planes defined by tuples of the two other principal axes (as shown in Figure 12), are shown in Figure 14. In this case, the rotation around ∣i∣ is rotation measured on the (j, k) plane with (i ≠ j ≠ k), as defined by a right-handed Cartesian coordinate system. The galactocentric radius is defined as ${r}_{{jk}}=\sqrt{{j}^{2}+{k}^{2}}$. In both figures, histograms in the right panels show the total velocity distribution (black), as well as the velocity distribution near the center of the galaxy (red). Stars in the center were selected based on their distance to the COM along each given axis using radii of 0.5, 0.33, and 0.3 kpc for the longest, intermediate, and shortest principal axis, respectively (following the axis ratios of the system).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The outer regions of Sgr are clearly affected by tidal forces. This can be observed as distortions in all of the rotation curves in the outer regions of the galaxy (see, for example, ${v}_{\mathrm{rot},x^{\prime} }$ in Figure 14, which shows positive rotation of the core up to distances ${r}_{y^{\prime} x^{\prime} }\sim 0.75\,\mathrm{kpc}$). Figure 13 shows the expansion of the galaxy along its longest principal axis, $x^{\prime}$, while the galaxy remains close to stable along $y^{\prime}$ and shows apparent contraction along $z^{\prime}$ for distances larger than ∼0.25 kpc. However, this apparent contraction seems to be the projection of the tidal tails' kinematics along $| z^{\prime} |$ (see plane x'–z' in Figure 26), rather than a real contraction of the core of the galaxy. We confirmed these findings, calculating the variation of the radial components of the velocity as a function of distance along each axis through finite differences. Similar results were obtained in all axes with expansion starting at different distances. We also found that the inner core of the galaxy shows almost no net expansion or contraction within an ellipsoid of axes (0.5, 0.33, 0.3) kpc. This ellipsoid is rotating mainly about its intermediate principal axis, $y^{\prime}$, at an average ${v}_{\mathrm{rot},y^{\prime} }2.46\pm 0.22\ \mathrm{km}\ {{\rm{s}}}^{-1}$. Residual rotation can also be observed about the shortest (${v}_{\mathrm{rot},z^{\prime} }2.11\pm 0.22\ \mathrm{km}\ {{\rm{s}}}^{-1}$) and, to a lesser extent, longest (${v}_{\mathrm{rot},x^{\prime} }1.38\pm 0.22\ \mathrm{km}\ {{\rm{s}}}^{-1}$) axes, indicating the presence of precession and nutation due to the MW's tidal torques. We aligned this inner ellipsoid to its angular velocity vector and measured a maximum rotation velocity of v_rot =4.13 ± 0.16 km s⁻¹ within this ellipsoid. It is important to notice that the orientation and strength of this vector vary with radius as shearing tidal forces change the shape and kinematics of the galaxy.

The 3D renderings of Sgr aligned to its principal axes of inertia on a galactocentric frame can be found in Appendix F. An interactive version of this rendering can be found at https://www.stsci.edu/~marel/hstpromo/Sagittarius_GC.html.

Our results allow us to discern between at least two components: a main body with a spheroid shape and the tidal tails.

Previous studies have shown that Sgr has an elongated body with some depth along the line of sight. Using red clump stars, Ibata et al. (1997) found that its core is consistent with a prolate spheroid with ratios 3:1:1 for its major, intermediate, and minor axes. More recently, Ferguson & Strigari (2020) used RR Lyrae stars from Gaia DR2 and the OGLE collaboration to find that the Sgr body is likely to have a triaxial geometry with ratios 1:0.76:0.43. Our results sit somewhere in between these two works, with a triaxial geometry (1:0.67:0.60) but closer to prolate than what Ferguson & Strigari (2020) suggested. In terms of depth, our derived FWHM of 2.48 kpc agrees with that from the OGLE collaboration (Hamanowicz et al. 2016), with an FWHM of 2.42 kpc (see Appendix C.1.1).

Our derived inclination of the Sgr main body also differs from the Ferguson & Strigari (2020) results. They found the longest principal axis of the galaxy to be almost parallel to the plane of the sky ($-{4.9}_{-18.8}^{+17.5\circ}$), whereas our results suggest an inclination of 43° ± 6° for the longest and 2° ± 7° for the intermediate axes in the inner regions of the galaxy, r < 1 kpc. It is important to notice that the overall inclination of the system depends on the considered radius. The inclination flattens as we consider larger radii due to the galaxy's elongation along its orbital path but remains as large as ∼35° up to distances of 2.75 kpc. This means that the major principal axis is inclined ∼30° with respect to the MW potential direction (see Figure 10). Stars closer to the MW center are subject to stronger tidal forces and thus get stripped earlier from the Sgr body while gaining angular velocity along their orbit. This mechanism is responsible for the observed "S" shape, a feature whose orientation is tightly constrained by the position of the galaxy relative to the Galactic center and the direction of motion. While tidal tails also occur in the N-body with spherical profiles, only ellipsoids forming certain angles (30° ≲ α ≲ 60°) will feature this characteristic "S" profile. In turn, models that well reproduce the Sgr stream properties show this shape regardless of the shape and internal kinematics of the progenitor. The fact that both the ML and N-body models made qualitatively consistent predictions about this feature strengthens the results obtained in the present work.

Another intriguing aspect is the bar observed in the galaxy. Such a bar was also reproduced by our rotating N-body models at the first pericenter passage on its orbit around the MW, after t = 0.46 Gyr from the start of the simulation. As described in detail by Łokas et al. (2014) and Gajda et al. (2017), the formation of a bar is a characteristic intermediate stage in the process of the transformation of disky dwarfs into dwarf spheroidals. This prolate shape is maintained and again enhanced right after the second pericenter passage, which occurs at 1.2 Gyr. The orientation predicted by the ML models is reached soon after, at around 1.22–1.24 Gyr. Arguably, a third passage would help to better reproduce the observed properties of Sgr, since it could naturally explain some of the stream's farthest features and would reduce the rotation signal of the remnant. However, the present model would result in a remnant that is too spherical and almost nonrotating if it survives the third pericenter, or, if it is less strongly bound, it will disrupt soon after the second pericenter. A model that survived a third pericenter while keeping some internal rotation would require a new selection of structural parameters of the dwarf, addition of the gas, star formation, etc., which would involve a computationally prohibitive parameter survey. We thus decided to stop the simulation after the second passage. In any case, we emphasize that in order to form a strongly elongated shape similar to the one seen in the Gaia data, the dwarf should transform into a bar at some stage of its evolution, and for this to occur, it is necessary that the progenitor has the form of a stellar disk.

It is worth mentioning that the 3D shape derived by our ML model (an elongated shape with tidal tails) falls in with those predicted by different N-body models whose main goal was to reproduce the properties of the core of the galaxy (see, for example, Łokas et al. 2010; Vasiliev & Belokurov 2020).

Several works have previously analyzed Sgr's structural and kinematic properties, looking for signs of the ongoing destruction of the galaxy. Aiming to reproduce the stream properties, Peñarrubia et al. (2010) proposed a scenario in which the Sgr progenitor is a late-type rotating disk galaxy inclined 20° with its orbital plane. Łokas et al. (2010) focused on reproducing the properties of the core of the galaxy from a similar starting point. Both works presaged a rotating core remnant for Sgr, which later works were unable to confirm. Of special relevance to this discussion are Peñarrubia et al. (2011) and Frinchaboy et al. (2012). Both used samples of giant stars with v_los, concluding that the Sgr kinematics were well reproduced by those of a pressure-supported system, i.e., nonrotating. Interestingly, a preliminary version of the data presented by Frinchaboy et al. (2012) was used in Łokas et al. (2010) as a comparison with the rotating remnant N-body galaxy. The Gaia mission has brought new possibilities of studying Sgr in detail; in the recent contribution of Vasiliev & Belokurov (2020), the authors did not find any significant rotation signal from their best-fitting N-body model. It is worth noting that some residual rotation (<4 km s⁻¹) was found by Frinchaboy et al. (2012), although they concluded that it could be due to other factors rather than intrinsic rotation from Sgr's progenitor.

Many aspects are contributing to the observed velocity field, making this a difficult issue. Only having access to the full phase-space information allows us to untangle the most relevant ones.

• 1.  
  The contribution of the COM motion (or bulk motion) to the observed velocity field in the sky. Here Sgr is moving along its projected optical major axis toward negative values of λ with a tangential velocity of 254.8 ± −2.8 km s⁻¹ projected on the sky. This tangential motion projects over the line of sight as a positive gradient in the same direction of movement, in this case 4.54 km s⁻¹ deg^–1, assuming a constant distance of D = 25.97 kpc.
• 2.  
  The Sgr geometry and inclination with respect to the celestial plane. Triaxial system Sgr has its longest principal axis inclined 43° with respect to the celestial plane. That translates to an average of ∼1.15 kpc distance difference between the areas closer to us and the farthest ones along the surveyed area in this work. This gradient in distance reflects on the projected tangential components of the velocity (PMs), with stars that are farther away apparently moving slower than closer ones. In the particular case of Sgr bulk PMs, the difference in distance projects as ∼16.7 km s⁻¹ in tangential velocities.
• 3.  
  The internal kinematics of Sgr. This includes both residual kinematics from the Sgr progenitor and perturbations of these from an external potential, such as that induced by tidal forces from the MW or the Large Magellanic Cloud.

In this paper, we have derived the full phase space for more than 1.2 × 10⁵ member stars, which allows us to simultaneously account for all of these effects and study the Sgr internal dynamics with unprecedented detail. We have shown that the Sgr core rotates in the plane of the sky at 2.9 ± 0.2 km s⁻¹ counterclockwise. It does so while also rotating perpendicular to the plane of the sky, with its east side moving toward the Sun (upon subtraction of the systemic motion of the galaxy).

It is clear that the internal dynamics of Sgr is linked to its orbital history and interaction with the MW. All tested N-body models show similar velocity gradients in the outer regions of the galaxy, caused by the tidal forces from the MW. However, we found substantial differences in the inner parts of Sgr. Our results show that despite being ripped apart by tidal forces, the galaxy conserves an inner region that is not expanding and is rotating at v_rot = 4.13 ± 0.16 km s⁻¹.

The spherical N-body models tested in this work failed to reproduce these characteristics, resulting in flat velocity profiles in the central regions of the galaxy (plots of these models are shown in Appendix E). This was expected, since spherical systems are not subject to torque forces. A flattened system, on the other hand, changes the part of the galaxy that is closer to the MW as it orbits around it, suffering torques that would induce velocity gradients. This results in velocity gradients along axes other than the original axis of rotation, ${v}_{z}^{\prime}$. As the system evolves and stretches along its orbit, the shape and kinematics of the system change in such way that most of the original rotation is conserved about the intermediate principal axis ${v}_{y}^{\prime}$, while some rotation is also transmitted to the other axes. This can be seen as a negative gradient in ${v}_{\mathrm{rot},x^{\prime} }$ along ${r}_{y^{\prime} z^{\prime} }$ or a positive gradient in ${v}_{\mathrm{rot},z^{\prime} }$ along ${r}_{x^{\prime} y^{\prime} }$ in Figure 14 that were well reproduced by the rotating model but not the spherical models.

The flattening of Sgr would imply some initial rotation of the system, which would also contribute to the observed gradients. However, the detected rotation velocity is still 1 km s⁻¹ below the predicted value of our best rotating N-body model candidate (3.5 ± 0.4 km s⁻¹). As commented in the previous section, a possible explanation for the excess of rotation in the rotating model would be that the simulated galaxy suffered just two pericenter passages. During each passage, the progenitor is stirred and perturbed, heating its kinematics and removing rotation. The Sgr is likely to have suffered more than two pericenter passages; only models with three or more passages are able to reproduce the stream's full extension (Dierickx & Loeb 2017). Furthermore, three main star formation events have been observed in the MW's star formation history at the approximate times when such passages occurred (Ruiz-Lara et al. 2020). However, even with just two passages, the rotating model does a good job reproducing the closer parts of the stream (see Figure 5 in Łokas et al. 2010). Future high-resolution N-body models aimed at reproducing both the stream and the core with three or more passages may provide some insight.

Our data allow us to derive the 3D vector for the internal angular momentum of the galaxy. For stars within ${r}_{3{\rm{D}}}=1\,{\rm{kpc}}$, this vector forms an angle θ = 18° ± 6° with respect to the orbital angular momentum of Sgr, i.e. the normal vector to the orbital plane of Sgr, meaning that Sgr is in an inclined prograde orbit around the MW. The slight inclination of the axis of rotation with respect to the orbital plane could aid in the explanation of the bifurcation of the stream (Belokurov et al. 2006; Koposov et al. 2012), as different parts of Sgr get stripped with slightly different angles as the dwarf moves along its orbit. This possibility was investigated by Peñarrubia et al. (2010), who found that a retrograde orbit, θ = −20°, was the best inclination that reproduced the observed bifurcation in the south arm. The model was later shown to be unsuccessful in reproducing the internal kinematics of Sgr (Peñarrubia et al. 2011). On the other hand, while our rotating model, prograde and initially inclined 10° with respect to its orbital plane, well reproduces Sgr's internal dynamics, it did not create two clearly distinct arms in the stream but rather a single wider one. Higher initial inclinations were also tested, although they showed a poorer match between the observed and simulated angular momentum vectors.

The 3D positions reveal still more differences between observations and retrograde models, such as the observed bar. Dwarfs in prograde orbits tend to form tidally induced bars after pericenter passages. This is not the case for retrograde systems, whose stellar content remains disky along the entire orbit (Łokas et al. 2015).

Very recently, Vasiliev & Belokurov (2020) presented N-body models that formed an elongated structure that resembles a bar despite not having initial rotation. The authors argued that the lower concentration used in their models allowed tidal forces to perturb the system to form a bar-like structure. While the morphology could be reproduced in such a way, the authors found that the resulting dwarf would be tidally disrupted in its totality. Our results, on the other hand, show that Sgr conserves a small inner core that is not expanding, supporting a scenario in which the inner regions of the galaxy might survive this last pericenter passage. This is also the case for the rotating model, which was the best at reproducing the configuration predicted by the ML models and whose inner core remains tidally bound at least until the next pericenter passage (see Figure 6 of Łokas et al. 2010). Therefore, we do not confirm the complete disruption of the Sgr after its recent pericenter passage.

Differences in the evolution of our N-body model from that of Vasiliev & Belokurov (2020) are to be expected given the differences in how models were compared to observations; models in Vasiliev & Belokurov (2020) were fitted directly to Sgr's observables projected on the sky, whereas in this work, a small set of models were compared to 6D predictions from ML. We also note the qualitative, "by-eye" nature of the model selection in both cases, which may weight the various observables differently in each case. The fact that in our model, Sgr's inner regions remain tidally bound after this pericenter passage might indicate that some adjustments should be made to the Vasiliev & Belokurov (2020) model's central initial conditions. But, given our own observational errors and model degeneracies, we do not make any strong claim about this particular aspect of Sgr evolution.

C.1.1. Errors in the Training Sample and Estimation of Real FWHM

C.1.2. Different Models

As shown in Appendix C.1.1, almost one-third of the observed dispersion along the line of sight (1.16 kpc) could be caused by observational errors. The ideal model will regularize the fitting to the training list, averaging several stars, avoiding overfitting the errors, and providing the most accurate distance possible. Hence, the dispersion of the residuals from the cross-validation is not indicative of how well the model performs when predicting distances, and we decided to use the Spearman correlation coefficient, the dispersion in the predicted distances in the central regions of the galaxy, and the existence of artificial gradients in the predicted distances with respect to the training sample.

The SVR with linear or polynomial kernels and the KNR models tend to produce better results using fewer features, while the SVR with Rbf kernel, MLP, and RFR benefit from including extra information, such correlations between features. This behavior is somewhat expected, as SVR and KNR models are simpler and may not obtain useful information from the correlation features, which show a periodic behavior in the sky along the covered area. In all cases, best results were obtained by fitting at least four features (columns).

The best model overall was SVR with the Rbf kernel fitting all features, although just marginally. This model is computationally expensive and was very closely followed by the KNR, which is the fastest of the seven tested models. We thus decided to use the KNR until less than 100 stars were rejected during an iteration in our code. Then we shifted to the SVR model with the Rbf kernel until convergence. Our final test results obtained with this method were the same within the errors as if only the KNR model was used; thus, we decided to keep the KNR results for their simplicity and faster performance.

The number of neighbors, k, and the best combination of independent variables, x , were determined on every iteration through the NCV test. Figure 19 shows three examples of results when using k = 1, 10, and 100 on the final training list. Low values of k produce slightly steeper predictions versus ground truth (closer to a one-to-one relation) but at the cost of having too high variance. Given our prior knowledge of the training list's errors, we favored configurations producing results with a central distance dispersion around 2.48 kpc measured for HB stars. In other words, we favored bias versus variance in order to prevent the artificial dispersion observed in the training list from propagating to our predictions. Notice that, as mentioned before, the error limitations in the training list prevented us from using the one-to-one relation as an indication of how well the model is performing. In fact, the KNR with the best one-to-one relation between D_pred and D_true, k = 1, was the worst model in any of the performed comparisons.

Zoom In Zoom Out Reset image size

Depending on the iteration of the membership selection, the best results were obtained using k between 10 and 30 and x = (x, y, μ_x , μ_y , ϖ). Larger values of k tend to produce overly flattened results without much improving the variance. The KNR that produced the best results in the final iteration was with k = 15, although values ranging from k = 9 to 32 provided consistent results. Results removing ϖ or μ_y were also reasonably consistent, indicating that some correlation between these two variables might exist.

In Figure 20, individual conditional expectation (ICE) curves show the dependence between the KNR predictions and each feature in x . The position on the sky, (x, y), and PM along the west direction, μ_x , have the largest impact on the prediction. On average, there is also a nonmonotonic dependence of the prediction on parallax, ϖ, which is not intuitive. This dependence may stem from the complex relationship of the random and systematic errors in Gaia parallaxes and PM. In any case, the dependence is weaker than for the first three features, and models not using ϖ or μ_y produced similar results to our standard model.

Zoom In Zoom Out Reset image size

The bright stars in the v_los training sample show relatively small errors (<2 km s⁻¹) compared to the observed dispersion. This mitigates the risk of overfitting and allows us to use the sum of the squared residuals as indicative of how well the model is performing. On the other hand, the stars in this training sample are sparsely distributed along the sky, which poses a challenge to the models that have to interpolate the gaps in a sensible way.

We tested models with nonlinear kernels, such as second- and third-degree polynomial or radial basis functions, with mixed results. While in some cases, nonlinear kernels provided better predictions in the training set, they tended to produce overly large or low velocities for extreme values of x in the main sample. This behavior can be controlled by means of regularization, $z^{\prime}$, but introduces a large dependence of the result on this and other metaparameters of the fitting. When using linear kernels, on the other hand, SVRs are highly insensitive to the use of larger $z^{\prime}$ and thus to overfitting. The choice of a linear kernel is thus justified in the outer parts of the galaxy by the simplicity of the model and the readability of its results, although one should be aware of its limitations in reproducing the nonlinear features observed in the training set.

To overcome this, we decided to combine our best-performing linear and nonlinear learners into an SR. In an SR, a model is trained over the prediction of several learners through cross-validation of their results. They generally outperform the best of the individual learners but are computationally expensive to run. Our SR model used a nonlinear artificial neural network, a MLP, to combine the predictions of six regressors, three nonlinear (extra trees, support vector machine with a radial basis function, and Gaussian process) and three linear (lasso, support vector machine with a linear kernel, and an elastic net). We chose the best configuration for each of the learners by maximizing ${R}^{2}={(1-\sum {({y}_{\mathrm{true}}-{y}_{\mathrm{pred}})}^{2}/\sum (({y}_{\mathrm{true}}-\left\langle {y}_{\mathrm{true}}\right\rangle )}^{2})$ in a grid covering different values of their fitting metaparameters through an NCV. The experiment showed that the model provides robust predictions against mid-to-small metaparameter variations. Only shifting all metaparameters toward a specific behavior, e.g., to have more regularization or large variations on the architecture of the final MLP metaregressor, had a nonnegligible impact on the final predictions. This was somewhat expected, as the metaregressor leverages the input of the six independent models by choosing the best possible prediction, hence making it more stable against spurious bad predictions.

Figure 21 shows an example for three different MLP architectures in the metaregressor. The use of more neurons and extra hidden layers can help to better predict velocities. However, adding more layers or neurons can result in an overfitted model, producing predictions with needlessly large variances. The selected best architecture for the MLP varies slightly depending on the combination of the input models' metaparameters, and we found that the best 25% of our models produced indistinguishable results within the errors. Our final model in the last iteration had an architecture of (5, 95, 85, 1) and uses a ReLU activation function and an L-BFGS solver. The ICE curves showing the dependence of the prediction on each feature in x are shown in Figure 22.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The resulting model was able to reproduce nonlinear features while showing stability against high-order oscillations on areas where no information is available. However, although the model produced some dispersion in the predicted velocities, cross-validation and final tests showed that the reproduced dispersion is likely to be driven by random differences with the actual velocities rather than an accurate prediction of naturally dispersed values. This manifests as the inability of the model to reproduce the extreme values of the velocity in the left panels of Figure 21, failing to achieve a one-to-one relation. This was the also the case for all individual regressors, regardless of their degree of complexity and nonlinearity. These findings were totally expected, as velocity stochasticity cannot be reproduced star by star by nonphysically motivated models such the ones used here. Furthermore, despite having much better errors than the distance training list, the line-of-sight velocity errors in PMs are still high star by star (∼0.12 mas yr⁻¹), causing random variance in the residuals. This does not pose a problem to our analysis, where our aim is to study the general velocity trends in the galaxy, averaging stars in Voronoi cells containing ∼1000 stars. Our NCV tests show that the SR does not produce any undesired spatial bias that could otherwise affect our results.

We have conscientiously selected, configured, and trained our models in order to minimize possible biases and errors in their predictions. The models presented in this work were able to make accurate predictions corroborated by our extensive cross-validation tests. The results obtained in the gaia_source table shown in this paper are fully compatible with those obtained separately for our training samples, where only one of the variables has been modeled (distance for the v_los catalog and v_los for the distance catalog). Furthermore, results obtained by imposing a constant D on the v_los training list or a constant v_los on the distance training list also provided compatible similar results, indicating that the ML models used are able to make good predictions based on the training list they used to learn. Yet the ML models are limited by the quality of both the target data, y_pred, and the features used to train the model and make predictions on new data, x . The observational errors affecting these quantities negatively affect the models' prediction abilities. We believe that, in terms of the internal dynamics and structure of Sgr, no meaningful star-by-star prediction can be made with our models given the current limitations of the data. Furthermore, training data are likely to be affected by systematic errors. A model that well reproduced a feature, y_true, in the training list affected by systematics would be propagating those systematics to its final predictions, y_pred.

We have performed extensive tests in order to find and analyze systematics in the values of y_true in our training lists, finding the data to be highly homogeneous and of similar quality. We have also propagated all sources of errors we are aware of to our final results and averaged these within Voronoi cells of ∼1000 stars, allowing statistically meaningful comparisons between cells. Lastly, we have inspected "by-eye" the predictions from our best candidate ML models, which all provided qualitatively consistent results between them, with our training lists and previous studies in Sgr. Thus, our decisions on the plausibility of our solutions are ultimately based on our previous knowledge of Sgr and comparisons to physically motivated N-body models as a final check.

Standard uncertainties, or statistical errors, accompany each magnitude in the Gaia DR2 catalog. These measure the precision rather than the accuracy and are known to be underestimated by a factor of 1.10 (Arenou et al. 2018; Lindegren et al. 2018). We propagated with Monte Carlo all statistical errors through all coordinate transformations performed. These include errors in the coordinates of the stars and the determination of the COM, σ(α*), σ(δ), $\sigma ({\alpha }_{0}^{* })$, and σ(δ₀), as well as errors in the PMs, σ(μ_α*), σ(μ_δ ), σ(μ_α*0), and σ(μ_δ0). Errors from external quantities such the distance to the system, σ(D₀), or its v_los, σ(v_los), are also considered and properly propagated.

The procedure for deriving distances and line-of-sight velocities is also performed in a Monte Carlo fashion, considering and propagating the errors from all independent and dependent variables. A total of 10⁶ realizations are performed, randomly sampling each independent and dependent variable from a normal distribution centered in its corresponding nominal value. The predicted values using the nominal values of the independent variables are assumed to be the nominal values of the final result. The quadrature addition of the standard deviation of all realizations plus the typical standard deviation from the cross-validation test is adopted as the total error.

Systematics are known to be affecting DR2 data, artificially shifting the astrometric zero-points. They are yet not well understood, but it appears that their effects can be differentiated into two components. The first one consists of large-scale systematics with a scale length of θ = 20°. With such a large scale length, the main effect these systematics will have in our data is to modify the bulk parallax and PMs of the galaxy. This will, in principle, have negligible effects on the detection of rotation.

We have studied any possible large-scale systematics effects on the bulk PMs and rotation using quasars. We found 278 AllWISE AGN sources (Secrest et al. 2015) that cross-match with Gaia DR2 sources within our field. We cleaned this sample of spurious sources by applying the procedure explained in Section 2, i.e., rejecting sources not compatible with the centroid of the distribution up to 3σ. The median of their PMs was added to the PMs of the stars in our sample in Table 1 and the subsequent calculations.

The small-scale systematics, on the other hand, introduce a variation in the zero-point of the parallax and PMs with a regular pattern of scales of ∼1° likely related to the Gaia scanning law. This banding pattern should average out for sufficiently large objects on the sky, such as Sgr. We have run tests to quantify how much these small-scale systematics may be affecting our measurements. Figure 23 shows the observed small variations in the astrometric zero-points caused by the small-scale systematics present in Gaia after removing most of the Sgr star. Trying to correct for these systematics could introduce an extra source of error. We therefore decided to calculate and propagate the variance introduced by these patterns in the final solution.

Zoom In Zoom Out Reset image size

Following the same prescription adopted by van der Marel et al. (2019), we create a mock systematic pattern with an rms of 35 μas yr⁻¹ (consistent with Gaia Collaboration et al. 2018b). We then shift the origin in sky coordinates (x, y) by random amounts between [−05, 05) and rotate the pattern between [−π, π). This is done separately for μ_α* and μ_δ 10⁶ times, recalculating the COM PM and the rotation signal about it in each iteration. From all realizations, we obtained a standard deviation of 0.0006 mas yr⁻¹ for both σ(μ_α⋆) and σ(μ_δ ) from the distribution of the COM PMs and 4 × 10⁻⁶ mas yr⁻¹ for σ(μ_r ) and σ(μ_ρ ). This supposes a variation of σ(v_W ) = 0.8 and σ(v_N ) = 0.7 km s⁻¹ and σ(v₃) =0.0005 and σ(v₂) = 0.001 km s⁻¹. These errors are included in all error intervals provided in this paper.