The last breath of the Sagittarius dSph
MMNRAS , 000–000 (0000) Preprint 5 June 2020 Compiled using MNRAS L A TEX style file v3.0
The last breath of the Sagittarius dSph
Eugene Vasiliev , (cid:63) , Vasily Belokurov Institute of Astronomy, Madingley road, Cambridge, CB3 0HA, UK Lebedev Physical Institute, Leninsky prospekt 53, Moscow, 119991, Russia
ABSTRACT
We use the astrometric and photometric data from
Gaia
Data Release 2 and line-of-sight velocities from various other surveys to study the 3d structure and kinematics ofthe Sagittarius dwarf galaxy. The combination of photometric and astrometric datamakes it possible to obtain a very clean separation of Sgr member stars from theMilky Way foreground; our final catalogue contains 2 . × candidate members withmagnitudes G <
18, more than half of them being red clump stars. We construct andanalyze maps of the mean proper motion and its dispersion over the region ∼ × ∼ ◦ from the orbitalvelocity and extending up to ∼ N -body simulations of a disrupting Sgr galaxy as it orbits the MilkyWay over the past 2.5 Gyr, which are tailored to reproduce the observed properties ofthe remnant (not the stream). The richness of available constraints means that onlya narrow range of parameters produce a final state consistent with observations. Thetotal mass of the remnant is ∼ × M (cid:12) , of which roughly a quarter resides in stars.The galaxy is significantly out of equilibrium, and even its central density is below thelimit required to withstand tidal forces. We conclude that the Sgr galaxy is terminallywounded, and will be fully disrupted over the next Gyr. Sagittarius dwarf galaxy (Sgr) is one of the closest and mostmassive satellites of the Milky Way (MW), however, itsstructure and properties are still not well known. By farthe most spectacular feature of this galaxy is the giant tidalstream covering a large fraction of the sky, indicating anongoing disruption of the satellite. The core of Sgr galaxyitself is located behind the Galactic bulge, and was discov-ered only relatively recently (Ibata et al. 1994, 1995). Due tothe small distance from the Galactic centre ( (cid:46)
20 kpc) and adistorted shape, it was immediately suspected to be tidallydisrupting, and indeed a few years later the full extent of thetidal stream was uncovered in the 2MASS survey (Majewskiet al. 2003). Since the stream spans a large range of galacto-centric radii and wraps around the Galaxy more than once,it has been studied extensively as a probe of the Galacticgravitational potential (e.g., Law & Majewski 2010). On theother hand, the dynamical state and properties of the Sgrremnant did not receive a comparable level of attention.Pe˜narrubia et al. (2010) developed a scenario in whichthe Sgr galaxy had an initially rapidly rotating stellar disc,which explained the bifurcation seen in the Sgr stream (Be-lokurov et al. 2006). However, subsequent observation of theline-of-sight velocity field across the remnant (Pe˜narrubia etal. 2011) did not agree with the predictions of that scenario.At the same time, (cid:32)Lokas et al. (2010) presented anothermodel for the Sgr remnant, which also started off as a ro- tating disc galaxy, but developed a strong bar induced bytidal torques from the MW, nearly eliminating all rotationin the remnant. This model was in a better agreement withthe then-available line-of-sight kinematics (Frinchaboy et al.2012), and remains the most recent model of the Sgr core.The original mass of the Sgr progenitor and the currentmass of the remnant have continued to be a subject of de-bate since the discovery of the dwarf. It became apparentsufficiently early on that a satellite as a large as Sgr comingas close to the Galactic disc could impart plenty of dam-age (see e.g., Ibata & Razoumov 1998; Bailin 2003). Whilethese earlier studies assumed a dwarf with a mass of orderof 10 M (cid:12) , the subsequent census of the stellar content ofthe satellite revealed that its total mass could be as high as10 M (cid:12) (see Niederste-Ostholt et al. 2010). This heavierSgr appeared a much more serious threat to the integrityof the Galaxy: now it could make the MW ring, seed spiralwaves (see Purcell et al. 2011) or even disrupt the stellardisc altogether (see Laporte et al. 2018).The second data release (DR2) of the Gaia space obser-vatory (Gaia Collaboration 2018a) dramatically expandedour knowledge of stellar kinematics in the MW and re-vealed that the Galactic disc is presently out of equilibrium,strongly implying a recent interaction with a massive per-turber, quite likely the Sgr dwarf (see e.g. Antoja et al. 2018;Carrillo et al. 2019; Laporte et al. 2019; Bland-Hawthorn etal. 2019). It remains to be established, however, whether c (cid:13) a r X i v : . [ a s t r o - ph . GA ] J un E. Vasiliev & V. Belokurov p C u m u l a t i v e nu m b e r o f s t a r s N ( < p ) astrometry only (primary sample)all criteria (primary sample)all criteria (all stars) Figure 1.
Distribution of membership probabilities according tovarious criteria: number of member stars with membership prob-ability not exceeding the value on the abscissa. Dotted green –primary sample (2 . × stars, ∼ × members), using astro-metric information only (parallaxes and PM); solid red – the sameprimary sample but using all available criteria (astrometry, CMDand location); dashed blue – entire sample (10 stars, 2 . × members). The use of multiple classification criteria results in avery sharp distinction between members and field stars: half ofmember stars in the primary sample have probability of member-ship p (memb) (cid:38) . p (memb) (cid:38) p (memb) (cid:38) these numerical models of the Sgr–MW interaction satisfyall of the observational constraints on the dwarf’s mass andits loss rate. While several studies used Gaia
DR2 data torefine the all-sky view of the Sgr stream (Antoja et al. 2020;Ibata et al. 2020; Ramos et al. 2020), no detailed analysisof the stellar kinematics in the Sgr remnant existed – anomission we aim to address in the present paper. This willpave the way to establishing the current dynamical state ofthe remnant, its bound mass and will help us constrain theSgr in-fall conditions.We first describe the procedure for selecting candidateSgr members from the entire
Gaia catalogue in Section 2,which simultaneously provides the kinematic maps of themean proper motion (PM) and its dispersion across thegalaxy. In Section 3 we augment these data with the line-of-sight kinematics derived from several existing datasets. InSection 4 we discuss the 3d structure of the Sgr remnantinferred from photometry, and estimate its stellar mass byexamining the distribution of stars in absolute magnitudes.We analyze the observed kinematic properties of the rem-nant in Section 5, and compare them with tailored N -bodymodels of a disrupting satellite in Section 6. Finally, Sec-tion 7 wraps up. We use a multi-stage procedure to select candidate memberstars and to determine their astrometric parameters. Theentire parent sample contains 1 . × stars with G < ◦ ≤ l ≤ ◦ , − ◦ ≤ b ≤ − ◦ and Declination − ◦ < δ < − ◦ , which have astrometric and colour information in Gaia . Wecompute the extinction-corrected G -band magnitudes ( G )and colours ( G BP , − G RP , ) as per Equation 1 in GaiaCollaboration (2018b). In the output catalogue, we retainall stars with 13 ≤ G ≤
18 ( ∼ ). This magnituderange includes blue horizontal branch (BHB), red clump(RC) and RR Lyrae stars in Sgr dSph. Note however thatthe uncertainties of these stars’ PM measurements are rel-atively large and they overlap strongly wit the MW fore-ground population in the color-magnitude diagram (CMD).On the other hand, Sgr red giant (RG) stars with G < ∼ . × stars, which is used in the initial classifica-tion and determination of the PM field of the dwarf. Thesestars must satisfy 13 ≤ G ≤ | (cid:36) | ≤ × (cid:15) (cid:36) and a num-ber of additional quality filters recommended by Lindegrenet al. (2018): astrometric excess noise < RUWE < . phot bp rp excess factor < . . bp rp .We build a multidimensional mixture model forthe twocomponents: Sgr (object) and MW (foreground or field)stars, writing the distribution function of c -th componentin terms of the following parameters: position on the skyplane ( α, δ ), parallax and PM ( (cid:36), µ α , µ δ ), position in the Gaia colour-magnitude diagram (CMD) G , G BP , − G RP , ,and – for a small fraction of stars – additional criteria suchas the line-of-sight velocity v LOS . The distribution functionof c -th component in each of these subspaces is denoted as p (loc)( c ) (location), p (ast)( c ) (astrometry), p (CMD)( c ) (photometry), p (vel)( c ) (line-of-sight velocity). These probability distributionsare normalized so that the integral of p ( ... )( c ) over its respectivesubspace is unity.In a given subset of stars (e.g., a specific region on thesky), the fraction of Sgr members η can be derived from thefollowing argument. Consider a star with observed param-eters X ≡ { X (loc) , X (ast) , X (CMD) , . . . } The likelihoods offinding such a star among Sgr members or among field starsare, respectively, L (memb) ,i = η p (memb) ( X i ) , L (field) ,i = (1 − η ) p (field) ( X i ) , (1)where p ( c ) ( X i ) is the product of distribution functions ineach subspace p (ast)( c ) ( X i ) × p (CMD)( c ) ( X i ) × . . . , and the prob-ability of this star to be a Sgr member is p (memb) ,i = L (memb) ,i L (memb) ,i + L (field) ,i . (2)On the other hand, the membership fraction of the entiresample is η = 1 N N (cid:88) i =1 p (memb) ,i . (3)It can be determined iteratively, by repeating the steps (1–3) until the estimate of η converges. On the other hand, thisiterative procedure may be also used to update the proba-bility distributions of both member and field stars. Namely,at each iteration, one recomputes the parameters of p ( ... )( c ) from the values of X i | Ni =1 , weighted by the current esti-mates of membership probability of each star. This approach MNRAS , 000–000 (0000) nternal kinematics of the Sagittarius dSph s u r f a c e d e n s i t y [ s t a r s / s q . d e g r ee ] Sgr members (major axis)Sgr members (minor axis)field stars 0 5 10 15 20distance [deg]020406080100120140 c u m u l a t i v e nu m b e r o f s t a r s × Figure 2.
Left panel: surface density profiles of Sgr members and field stars as functions of distance. For Sgr members, we plot separatelythe profiles along the major (solid blue) and the minor (dot-dashed green) axes as functions of the distance from the dwarf’s centre. Forfield stars (dashed red), the abscissa is the offset in Galactic latitude b from the region boundary at b = −
6. The density of Sgr members(except for the central bump, corresponding to the globular cluster M 54) is well described by a King profile with axis ratio ∼ ◦ , and King parameter W = 4, which has a tidal radius of ∼ ◦ . Right panel: cumulative number of Sgr members as a function of distance. Solid blue and dashed magenta lines show the profiles alongthe major axis in the direction of the trailing and leadin arm, respectively. The latter one is truncated at ∼ ◦ since the footprint of ourcatalogue ends there, but it already shows signs of incompleteness at smaller distances. Dot-dashed green line show the profile along theminor axis (symmetrized and divided by two) as a function of distance multiplied by the axis ratio 2.8. is known as the expectation/maximization (EM) algorithm(e.g., Dempster et al. 1977).We follow this strategy, splitting it into several stages,in which the distribution functions in different subspacesare updated one by one. At the first stage, we determinethe distribution functions of member and field stars in the3-dimensional astrometric subspace (parallax (cid:36) and PM µ α , µ δ ). These are represented as a 5-component Gaussianmixture model, with the narrowest component being the Sgrstars, and the remaining ones represent the foreground pop-ulation. We use the Extreme Deconvolution method (Bovyet al. 2011), which takes into account the observational er-rors of each star and constructs the intrinsic distributionfunction, which is broadened by measurement errors beforeevaluating the likelihoods (1). We increase the parallax andPM errors quoted in the Gaia catalogue by 10% to compen-sate a slight underestimate of formal errors, as discussed byLindegren et al. (2018). Since the PM distributions of boththe Sgr and the MW stars vary considerably across the en-tire region of the sky, we perform this analysis separately in8 partially overlapping macro-regions on the sky (4 strips in b and 2 strips in l ). Already at this stage it became clear thatone needed to take into account the spatial gradients in thePM of the Sgr stars, therefore we subtracted a global lineartrend from the PM of all stars before running the mixturemodel on the shifted PMs µ (cid:48) α , µ (cid:48) δ defined as µ (cid:48) α ≡ µ α − . α − α ) + 0 . δ − δ ) ,µ (cid:48) δ ≡ µ δ + 0 . α − α ) + 0 . δ − δ ) , (4)where α = 283 . ◦ , δ = − . ◦ are the coordinates ofthe Sgr centre.The first step (based on astrometry alone) already al-lows one to build a representative sample of Sgr members(Figure 1, dotted green line). In the next step, this astro-metric classification is used to determine the spatial distribu- tion and construct colour–magnitude diagrams (CMDs) ofboth populations in the two-dimensional space of extinction-corrected apparent magnitudes G vs. G BP , − G RP , . Thespatial variation of the stellar density of the MW field starsis assumed to depend on the Galactic latitude b only, whilethe density of the Sgr member stars is assumed to be afunction of ellipsoidal radius ˜ R ≡ (cid:112) X + ( Y /q ) , where X and Y are the coordinates along the major and minor axes,respectively, and q is the flattening parameter. Both den-sity profiles are represented as free-form cubic splines with ∼
10 nodes each. We account for the fact that the foot-print of our catalogue is truncated at ∼ ◦ from the Sgrcentre in the direction of the leading arm by doubling thecontribution of all stars at distances > ◦ in the trailingarm to the density profile of members. The CMD distribu-tion functions are represented as 2d histograms in the range13 < G < , − < G BP, − G RP, <
4; we use 0.05-mag bins and additionally smooth the histograms with a0.1-mag Gaussian kernel. As the photometric errors are rel-atively small and depend mainly on the G -band magnitude,the CMDs are constructed for the error-broadened, not in-trinsic distributions (unlike the astrometric distributions).We follow the iterative EM procedure outlined above. Ateach iteration, the CMD histograms, the values of splinefunctions and the parameter q are recomputed, using thecurrent membership probability estimates p (memb) ,i of eachstar as weight factors. Then these probability estimates areupdated using the improved distribution functions of bothpopulations (Equations 1, 2), and finally the overall fractionof the Sgr members η is updated in Equation 3. Initially, wedetermine density profiles and global CMD for both pop-ulations in the entire region, and then proceed by build-ing more localized CMDs of field stars in separate macro-regions, while keeping the overall density profiles and the MNRAS000
4; we use 0.05-mag bins and additionally smooth the histograms with a0.1-mag Gaussian kernel. As the photometric errors are rel-atively small and depend mainly on the G -band magnitude,the CMDs are constructed for the error-broadened, not in-trinsic distributions (unlike the astrometric distributions).We follow the iterative EM procedure outlined above. Ateach iteration, the CMD histograms, the values of splinefunctions and the parameter q are recomputed, using thecurrent membership probability estimates p (memb) ,i of eachstar as weight factors. Then these probability estimates areupdated using the improved distribution functions of bothpopulations (Equations 1, 2), and finally the overall fractionof the Sgr members η is updated in Equation 3. Initially, wedetermine density profiles and global CMD for both pop-ulations in the entire region, and then proceed by build-ing more localized CMDs of field stars in separate macro-regions, while keeping the overall density profiles and the MNRAS000 , 000–000 (0000)
E. Vasiliev & V. Belokurov BP − RP G Gaia RVSApogeePenarrubia+ 2011Frinchaboy+ 2012 10 5 0 5 µ α µ δ Figure 3.
Left panel: extinction-corrected CMDs of Sgr members (solid contour lines) and MW field stars (dashed), as inferred fromthe combination of all classification criteria. Points show the stars with measured line-of-sight velocities from various datasets. Contoursare separated by 0.5 dex. The bright part of the Sgr red giant branch has virtually no overlap with the MW stars, which was used inprevious studies to select Sgr members without relying on astrometric information. However, for the majority of Sgr stars in the finalcatalogue, photometry alone would not be sufficient to determine membership, although it does help to sharpen the classification.
Right panel: distribution of Sgr members (solid contours) and MW field stars (dashed) in the PM space; contours separated by 0.5 dex.
Sgr CMD fixed. The resulting density profiles, CMDs andPM distributions are shown in Figures 2, 3.After constructing the CMD distribution functions, wereturn to estimating the parameters of the astrometric dis-tribution function p (ast)(memb) of Sgr members, but now deter-mine the mean PM and its dispersion tensor in smaller spa-tial regions. We divide the entire region into 200 polygo-nal bins containing roughly equal numbers of member stars,using the Voronoi binning scheme of Cappellari & Copin(2003) with some manual postprocessing. For all stars ineach region, we again run the EM procedure, at each itera-tion recomputing the 5 parameters of the Gaussian distribu-tion p (ast)(memb) : mean µ α , µ δ and the components of symmetriccovariance matrix. At this stage, we use the CMD distribu-tions of both populations determined previously, but not thedistance distributions (since these are normalized to unityin the entire region, not in each bin). We also keep fixed theparameters of the field distribution p (ast)(field) (recall that theyare also spatially-dependent, but vary on a larger scale thanthe bin sizes). We use only stars in the primary sample todetermine the parameters of p (ast)(memb) , but then compute themembership probabilities for all stars. This avoids a possibleinflation of the PM dispersion by fainter stars with lower PMprecisions, and more importantly, by stars with unreliableastrometry, which did not pass the quality filters. We esti-mated the uncertainties on the mean PM (which are domi-nated by the spatially correlated systematic errors in Gaia astrometry), using the method detailed in Vasiliev (2019b);these are (cid:46) .
05 mas yr − . The statistical uncertainties onthe PM dispersions are (cid:46) .
01 mas yr − . The resulting kine- matic maps (mean PM and its dispersion) are discussed inSection 5, and we provide the derived values in Table A1.Finally, we again run the EM procedure on the entiresample, keeping fixed all distributions except the densityprofile, to obtain the membership probability for each star(2). Figure 1 demonstrates that the combination of all se-lection criteria produces a very sharp distinction betweenmembers and field stars (almost 90% of candidate memberstars among the primary sample have membership probabil-ity exceeding 90%, while this fraction drops to 75% for theentire catalogue). The right panel of Figure 2 shows that thethe number of Sgr members in the leading arm (closer to theGalactic plane) is ∼ −
15% lower than in the trailing armwithin the same distance from its centre. This indicates thatour catalogue is somewhat incomplete at low Galactic lat-itudes, but it remains very pure: the sheer excess of fieldstars means that only very few actual Sgr members havea high enough likelihood ratio to be classified as such. Thisconstrasts with a more traditional selection procedure basedon fixed boxes in the CMD and PM spaces, which becomesmore contaminated as the density of field stars increases.The entire catalogue of candidate members is available in theelectronic form at https://zenodo.org/record/3874830 . Some of the brightest Sgr members have line-of-sight veloc-ity measurements from the
Gaia
RVS instrument ( ∼ MNRAS , 000–000 (0000) nternal kinematics of the Sagittarius dSph consider three complementary catalogues. Pe˜narrubia et al.(2011) observed ∼ ◦ from theSgr centre along the major axis and up to 2 ◦ along the minoraxis; almost 90% of them are actual Sgr members. Frinch-aboy et al. (2012) observed ∼ − ◦ to 12 ◦ along the majoraxis, − ◦ to 5 ◦ along the minor axis, and a few diagonallysituated fields. Roughly a half of these stars are Sgr mem-bers. Finally, several fields of the APOGEE spectroscopicsurvey are within the footprint of Sgr, and contain ∼ Gaia
RVS stars, which are generally brighter and more scat-tered across the region). We find that the v los measurementsare largely consistent between independent datasets withinerror bars. The differences between velocities of individualstars are of order a few km s − , and the systematic offsetsbetween entire samples are at a level of 1 − − , sig-nificantly smaller than the velocity dispersion of Sgr or thegradient of the mean v los across the galaxy. We thus com-bine the information from all available sources, for a totalof 3300 member stars with v los measurements.We group this spectroscopic sample into 36 Voronoibins, with the majority (30) being located within ∼ ◦ fromthe Sgr centre, and remaining ones (mostly from the Frinch-aboy et al. 2012 sample) spread along the trailing arm up to12 ◦ from the centre. We compute the mean v los and its dis-persion in each bin, taking into account the measurementuncertainties. Since the values of v los for individual starswere taken into account together with other astrometric andphotometric properties when determining the membershipprobability, we do not need to impose further filters on thesample. The statistical uncertainties on the mean v los are atthe level 1 − − , depending on the number of stars inbins. We provide the derived values in Table A2 and discussthe kinematic maps in Section 5. It is apparent that the Sgr remnant is stretched along its or-bit in projection, however, its 3d structure is less well known.Ibata et al. (1997), based on the then-available photomet-ric data, examined the magnitude distribution of Sgr starsaround the red clump (RC), which serves as an absolutephotometric reference point. They found no significant gra-dient of the RC magnitude (equivalently, variation of thedistance) across the face of the galaxy, and adopted the dis-tance D = 25 ± ∼ .
04 mag,which corresponds to the upper limit on the thickness of ∼ . , comparable to the width along the projectedminor axis. Hence they concluded that the Sgr remnang hasa strongly prolate cigar-like shape with axis ratios ∼ D (cid:39) . σ ∼ ∼
800 RR Lyrae stars from the
Gaia catalogue of variable stars, to study the 3d shape and orien-tation of the Sgr core, assuming it to be a triaxial ellipsoidwith a Gaussian density profile. They concluded that theintrinsic shape of the Sgr core is indeed triaxial, with axisratios 1 : 0 .
76 : 0 .
43, with the major axis nearly perpendic-ular to the line of sight, and the intermediate axis parallelto it. Their estimate of the distance to the Sgr centre isslightly smaller, 26.4 kpc, and the length (dispersion of theGaussian) of the major axis is 1.7 kpc. We adopt a distance D = 26 . Gaia photometry.For the RC stars, which constitute roughly half of our finalsample, we determined the peak and width of their distribu-tion as follows. We selected all stars with 17 . < G < . < G BP , − G RP , and membership probability above0.8. Then for each of the 200 Voronoi bins, we constructa two-component truncated Gaussian mixture model for thedistribution of stars in G , taking into account the sharpboundaries of the selection box (Figure 4, left panel). Thedispersion of the narrower Gaussian component ranges from ∼ .
06 in the centre to ∼ .
15 in the trailing arm, corre-sponding to the thickness along the line of sight ∼ . − . (cid:38) ◦ from the centre, the meandistance to stars starts to decrease. This apparently corre-sponds to the transition from the cigar-shaped remnant tothe trailing arm of the stream, which has a minimum helio-centric distance ∼
20 kpc at ∼ ◦ from the centre. We cau-tion that the variation in apparent magnitudes is fairly small In this section, we refer to the dispersion σ of the Gaussiandistribution as the thickness; the full width at half maximum is,as usual, 2.35 σ .MNRAS000
20 kpc at ∼ ◦ from the centre. We cau-tion that the variation in apparent magnitudes is fairly small In this section, we refer to the dispersion σ of the Gaussiandistribution as the thickness; the full width at half maximum is,as usual, 2.35 σ .MNRAS000 , 000–000 (0000) E. Vasiliev & V. Belokurov G d N / d G all starsRC starsother 10 -5 -4 -3 -2 N( Figure 4. Left panel: distribution of extinction-corrected apparent magnitudes of ∼ ◦ located 4 ◦ from the centre ( α = 288 . ◦ , δ = − . ◦ ). Solid blue line shows the distribution of all stars, which is fit by two truncatedGaussians in the range 17 . < G < 18, shown by green dashed and red dot-dashed lines. We assume that the narrower Gaussiancorresponds to RC stars, and use its mean and dispersion to estimate the photometric distance and thickness (this analysis is performedseparately for each of the 200 Voronoi bins, containing ∼ Centre panel: CMD of Sgr (gray) compared to those of two globular clusters with comparable metallicity [Fe/H] ∼ − . ± . 2: 47 Tuc(NGC 104, green) and NGC 6723 (red). Magnitudes and colours are extinction-corrected and translated to the absolute magnitudes. Right panel: cumulative number of stars brighter than the given absolute magnitude, normalized by the mass of the stellar system. Themasses for the two globular clusters are determined by Baumgardt et al. (2019) from dynamical modelling, while the mass of the Sgrremnant is chosen to match the other two curves, and is estimated to be in the range (1 − . × M (cid:12) . ( (cid:46) . − . 9; for comparison, wechose two clusters with large enough number of stars andlow reddening – NGC 104 (47 Tuc) and NGC 6723. Theseclusters have been extensively studied, and in particular,Baumgardt et al. (2019) fitted N -body dynamical models tothe cluster kinematics and photometry, and estimated theirtotal masses. We may infer the stellar (but not dynamical)mass of Sgr remnant by comparing the distribution of itsstars in absolute magnitudes to these clusters, after normal-izing the star counts by total masses. Figure 4, right panel,shows that a good match of the cumulative number of starsas a function of absolute magnitude is obtained for the stel-lar mass of Sgr around 10 M (cid:12) . We stress that this esti-mate does not assume that Sgr is in dynamical equilibrium,nor that its stellar mass is equal to the dynamical mass.These assumptions are only made for the globular clusters,and in addition we assume that the stellar mass functionsof Sgr and the clusters are similar. Our approach neglectsthe fact that stars in Sgr are generally younger than in theglobular clusters chosen for comparison. We did not calcu- late the total luminosity of the Sgr remnant, but if its stel-lar mass-to-light ratio is also similar to that of the clusters( M/L V (cid:39) . . × L (cid:12) . For comparison, Niederste-Ostholt et al. (2010)estimated it to be (0 . ± . × L (cid:12) , and Majewski etal. (2003) quote a value twice as smaller. Their estimate isbased on a different approach – extrapolating the surfacebrightness profile after subtracting the MW foreground. The Sgr centre, which coincides with the globular clusterM 54, has coordinates α = 283 . ◦ , δ = − . ◦ . Weadopt the distance D = 26 . v los , = 142 km s − , and PM µ α, = − . − , µ δ, = − . 35 mas yr − . We adopt the following values for the Solarposition and velocity in the Galactocentric rest frame: X = − . V { X,Y,Z } = { . , . , . } km s − (Astropycollaboration 2018); the corresponding position and veloc-ity of the Sgr centre are { X, Y, Z } = { . , . , − . } kpc, V { X,Y,Z } = { . , − . , . } km s − .Before presenting the kinematic maps, we need to de-fine the coordinates and quantities to plot. The analysis de-scribed in Section 2 used the PM in equatorial coordinates α, δ ; however, the mean PM and their dispersion tensor canbe equivalently transformed into any other celestial coordi-nates of choice χ, ξ , using standard expressions for sphericalgeometry. It is natural to assign χ = ξ = 0 to the Sgr cen-tre, but the orientation of the axes remains a free parameter.We now argue that it makes sense to align one of the coor-dinate axes (say, χ ) with the direction of the mean PM of MNRAS , 000–000 (0000) nternal kinematics of the Sagittarius dSph the Sgr core on the sky plane, so that µ χ, = (cid:113) µ α, + µ δ, and µ ξ, = 0. This direction is different from the major axisof Sgr (which itself very nearly coincides with its orbit onthe sky plane), because the observed PM has a contributionfrom the solar motion. However, it is in these coordinatesthat the perspective effects have the most straightforwardmanifestation.Consider a situation when an isolated, non-rotatinggalaxy with isotropic velocity dispersion σ is located at adistance D and moves with a velocity v relative to theobserver, directed parallel to the χ axis. The mean PM ofstars at a distance D is µ χ = v /D and µ ξ = 0, and its dis-persion is σ/D for both components. Since the galaxy hasa finite thickness, one needs to integrate along the line ofsight to obtain the PM dispersion σ µ . Assuming for defi-niteness that the density profile along the line of sight is aGaussian with a dispersion h (cid:28) D , it is easy to show thatin the first approximation, the average PM is µ χ, = v /D , µ ξ, = 0, and its dispersion is σ χ = (cid:112) σ + ( µ χ, h ) /D , σ ξ = σ/D . In other words, the PM dispersion along thedirection of motion is broadened by the spread in distances,while in the perpendicular direction it remains the same asif the galaxy had zero thickness. For Sgr, v (cid:39) 380 km s − , σ (cid:39) 13 km s − , h (cid:39) D = 26 . σ χ have comparable magnitudes. There-fore, the inflation of PM dispersion along the direction ofmotion due to perspective effects is very significant, and thealignment of the χ axis along the PM vector roughly diag-onalizes the PM dispersion tensor, justifying our choice ofthe coordinate system.Continuing with our toy example, if a galaxy movingas a solid body (with a spatially uniform mean velocity v )subtends a finite region on the sky, the observed PM and v los field will not be constant due to perspective distortions.Let v los , be the component of velocity along the line ofsight passing through the galaxy centre ( χ = ξ = 0), and v tan , ≡ µ χ, D be the velocity component parallel to the χ axis (the third component is µ ξ, = 0 by construction).Consider a star at a distance D = D (1 + ζ ), with thedimensionless parameter ζ (cid:28) χ, ξ (cid:28) 1, moving with a total velocity v + u . Its relativevelocity with respect to the galaxy centre u has components u χ , u ξ , u ζ in three perpendicular directions. To a first orderin χ, ξ and ζ , the observed PM and v los of this star are µ χ = µ χ, − v los , /D χ − µ χ, ζ + u χ /D,µ ξ = µ ξ, − v los , /D ξ − µ ξ, ζ + u ξ /D,v los = v los , + µ χ, D χ + µ ξ, D ξ + u ζ . (5)The terms proportional to χ, ξ are caused by perspec-tive effects: in the direction different from the galaxy centre,the line-of-sight velocity has a contribution from the centre-of-mass PM, and vice versa. The magnitude of these correc-tions is quite significant and can be larger than the relativevelocity components, but it involves known quantities, andcan be subtracted to obtain the ”perspective-corrected” PMfield: µ (cid:48) χ ≡ µ χ + v los , /D χ = µ χ, − µ χ, ζ + u χ /D,µ (cid:48) ξ ≡ µ ξ + v los , /D ξ = µ ξ, − µ ξ, ζ + u ξ /D. (6)However, the term proportional to ζ also has an equally significant contribution, and cannot be corrected since thedistance offset ζ is unknown. With our choice of µ ξ, = 0,this term is zero in the second row, so µ (cid:48) ξ does contain onlythe actual relative velocity field (still scaled by the unknowndistance).The line-of-sight velocity can also be corrected for theperspective effects, and does not contain any terms propor-tional to the unknown distance offset ζ . Obviosly, this ispossible only if the tangential component of the total ve-locity of the object’s centre of mass relative to the observer v tan , = µ χ, D is known. When the object’s PM is unavail-able, people often use a “partially corrected” quantity – so-called Galactic Standard of Rest (GSR) velocity v los , GSR . Itis defined as the velocity that would be measured by an ob-server residing at the Solar position, but with zero velocitywithin the MW: v los , GSR ≡ v los + v (cid:12) · n , (7)where v (cid:12) is the 3d solar velocity in the MW rest frame, and n ( χ, ξ ) is the unit vector in the direction of the given point { χ, ξ } on the celestial sphere. In other words, this correctioninvolves only the solar velocity, but not the total centre-of-mass velocity of the object, and hence does not get rid ofall perspective effects. By definition, v los , GSR measures thevelocity component in the direction of the observer , and thisdirection varies across the object. Thus, even if all stars ina galaxy were moving with the same 3d velocity, the ob-served v los , GSR would measure projections of this velocityonto different lines of sight, and hence would not be spa-tially uniform. Conversely, a constant v los , GSR field does notimply the absense of internal rotation. Therefore, v los , GSR does not have any advantages compared to v los for describ-ing the internal kinematics (in fact, a possible disadvantageis that the transformation between v los and v los , GSR obvi-ously depends on the adopted spatial velocity of the Sun,which may differ between studies). However, it turns outthat the mean value of v los , GSR varies only mildly across theSgr remnant, unlike either the heliocentric v los or the in-ternal velocity u ζ , which is nothing more than a fortuitouscoincidence. For this only reason, we will be plotting v los , GSR to highlight the small differences between observations andmodels, which otherwise would have been swamped by thestrong gradients caused by perspective effects. Figure 5 summarizes all observational data on the kinemat-ics of the Sgr core. Top row shows the mean perspective-corrected PM µ (cid:48) χ , µ (cid:48) ξ (panels A, B) and v los , GSR (panel C).As discussed above, µ (cid:48) ξ contains only the internal velocitycomponent u ξ /D ; it is reasonably flat over the main bodyof the dwarf and shows a steady gradient with χ at theedges of the remnant where the stars move away from thecentre into leading and trailing arms. The other component, µ (cid:48) χ , has contributions both from the internal velocity com-ponent u χ /D and from the distance gradient − µ χ, ζ . Thesharp increase in µ (cid:48) χ happening in the trailing arm ∼ ◦ from the Sgr centre indicates the sudden drop in the meandistance to stars, not the change in their internal velocity.This is corroborated by a similar feature in the photomet-ric distance map (panel J). Physically, this corresponds tothe transition between the Sgr remnant itself, which is tilted MNRAS000 Figure 5 summarizes all observational data on the kinemat-ics of the Sgr core. Top row shows the mean perspective-corrected PM µ (cid:48) χ , µ (cid:48) ξ (panels A, B) and v los , GSR (panel C).As discussed above, µ (cid:48) ξ contains only the internal velocitycomponent u ξ /D ; it is reasonably flat over the main bodyof the dwarf and shows a steady gradient with χ at theedges of the remnant where the stars move away from thecentre into leading and trailing arms. The other component, µ (cid:48) χ , has contributions both from the internal velocity com-ponent u χ /D and from the distance gradient − µ χ, ζ . Thesharp increase in µ (cid:48) χ happening in the trailing arm ∼ ◦ from the Sgr centre indicates the sudden drop in the meandistance to stars, not the change in their internal velocity.This is corroborated by a similar feature in the photomet-ric distance map (panel J). Physically, this corresponds tothe transition between the Sgr remnant itself, which is tilted MNRAS000 , 000–000 (0000) E. Vasiliev & V. Belokurov with respect to the line of sight and its own orbit, and theunbound tail, which is roughly parallel to the orbit, but liesat a larger distance. The heliocentric distance to the Sgrorbit decreases towards the trailing arm, but the mean dis-tance to stars in the remnant increases until the transitionzone.Artifacts from the Gaia scanning law manifest them-selves as systematically offset mean PM in spatial regions ∼ . ◦ across, most notably as a blue scar in the middle ofthe top left panel. The magnitude of these systematic errorsis (cid:46) . − , and does not obscure the real featuresseen in the data.The middle row displays the PM and v los dispersions.Panels D and E show the two components of PM disper-sion tensor in these coordinates, confirming our expectationsdiscussed above. The perpendicular component σ ξ is nearlyconstant (0 . − . 14 mas yr − ) across the field of view ,while the dispersion parallel to the direction of motion σ χ ranges from 0.15 to more than 0.35 mas yr − . We may usethe above toy example to estimate the “kinematic thick-ness” of the Sgr remnant from the difference between thetwo components of PM dispersion (assuming that the veloc-ity dispersion is isotropic, which, as we shall see, is not quitetrue): h kin ≡ D (cid:113) σ χ − σ ξ (cid:14) µ χ, . (8)This quantity is plotted on the panel G, and resembles qual-itatively the photometric thickness map (panel H), withlower values ∼ v los , GSR map obtained by interpolating among nearest 100stars at each location. There is a mild gradient of v los , GSR parallel to the major axis of Sgr, but as discussed above,this quantity does not have a straightforward physical inter-pretation by itself. Its dispersion, however, is a real phys-ically relevant quantity, and is remarkably constant acrossthe galaxy (panel F). σ los lies in the range 12 − 14 km s − ,and is slightly lower in the centre (Majewski et al. 2013).Remarkably, the PM dispersion σ ξ translated into km s − is very similar to σ los . However, one cannot conclusively in-terpret this as a sign of isotropy, since the third velocitydispersion component contributes only a fraction of the PMdispersion σ χ , and hence cannot be measured directly.However informative these kinematic maps are, they arestill not sufficient to reconstruct the internal velocity field u within the galaxy. Two of its components, u ξ and u ζ , canbe read off the panels B and C ( µ (cid:48) ξ and v los , GSR ); however,the PM component µ (cid:48) χ contains entangled information aboutboth the component of velocity u χ parallel to the apparentdirection of motion and the mean distance to stars, which isnot known to a sufficient accuracy to be subtracted. There-fore, a proper dynamical model is needed to interpret theobservations. N -BODY MODELS6.1 General considerations It is clear that the Sgr remnant is a heavily perturbed stellarsystem, and modelling it within the steady-state approxi-mation would be inadequate. Instead, we explore evolution-ary models of a disrupting satellite around the Galaxy, con-straining them to have the present-day position and velocityof the Sgr remnant.We set up an equilibrium model for the Sgr galaxy,as described below, using the Agama framework (Vasiliev2019a). We then place it roughly in the apocentre of its orbit ∼ M (cid:12) , according to Jiang & Binney 2000 orGibbons et al. 2017) than our adopted range of initial massat the time 2 Gyr ago ( ∼ M (cid:12) ), but the entire evolu-tion is outside the scope of the present study; we are onlyconcerned with the present-day state of the Sgr remnant.We evolve the Sgr galaxy under its own self-gravity plusthe static external tidal field of the Galaxy, assuming thatthe latter is fixed and not perturbed. We ignore the effectof dynamical friction and the response of the MW to thegravitational tug from its largest satellite – Large Magel-lanic Cloud (LMC). For our range of initial masses, dy-namical friction would change the orbit parameters by lessthan 10% per orbital period. Likewise, the reflex motion andthe distortion of the MW halo introduced by LMC (see e.g.Garavito-Camargo et al. 2019; Erkal et al. 2019) are impor-tant for modelling the Sgr stream (e.g. Vera-Ciro & Helmi2013; G´omez et al. 2015), but this is not the goal of thepresent study. Instead, we adopt a reasonably realistic MWmodel, with parameters drawn from an ensemble of MonteCarlo samples from McMillan (2017), but with a less massivespherical halo than in their best-fit model; the circular veloc-ity ranges from 225 to 185 km s − between 15 and 60 kpc. Inthis potential, the trailing arm of the Sgr stream aligns wellwith the observations, but the leading arm plunges back intothe Galactic disc too early; it is known that a single sphericalpotential cannot simultaneously fit both arms of the stream(see e.g. Helmi 2004; Johnston et al. 2005; Law & Majewski2010). However, it provides a good fit to the Sgr remnant.We run the simulations with the N -body code gyr-falcON (Dehnen 2000), which is included in the Nemo framework (Teuben 1995), with the external potential ofthe MW provided by the Agama plugin for Nemo . Thenumber of particles is a few × , softening length is (cid:15) =0 . 05 kpc (Plummer equivalent is 0.035 kpc), and the maxi-mum timestep is ∼ We explore many variants for the initial structure of theSgr progenitor. Single-component models were found to be MNRAS , 000–000 (0000) nternal kinematics of the Sagittarius dSph 15 10 5 0 55051015 [mas/yr] A [mas/yr] B v los , GSR [km/s] C [mas/yr] D [mas/yr] E - - los [km/s] F Kinematic thickness [kpc] G Photometric thickness [kpc] H Distance [kpc] J Figure 5. Kinematic maps of the Sgr remnant. Coordinates are aligned with the apparent (non-reflex-corrected) motion of the objecton the sky (the motion is in the direction of increasing χ ); the true velocity of the Sgr centre (blue arrow in panel G) points towards theGalactic plane, which lies slightly beyond the top right corner. A grid of equatorial coordinates is shown in Panel F. Most panels alsodisplay the surface density, with contours logarithmically spaced by 0.4 dex (i.e., one magnitude).Panels A and B show the perspective-corrected mean PM µ (cid:48) χ , µ (cid:48) ξ (Equation 6); panels D and E – the PM dispersions σ χ , σ ξ . The first ofthese values is inflated due to a non-negligible thickness of the galaxy, and the difference between the two dispersions (panel G) can beused to estimate the “kinematic thickness” h kin (Equation 8). It increases from ∼ ∼ v los , GSR (Equation 7) and its dispersion. Panel J shows the photometric distance estimate, which correlates with the featuresseen in the mean PM µ (cid:48) χ parallel to the direction of motion due to perspective effects. Colour scales for v los , GSR and σ los match thoseof PM for an assumed distance D = 26 . unable to match all observational constraints (see Niederste-Ostholt et al. 2012; Gibbons et al. 2017), so we concentrateon two-component models with more centrally concentratedstellar distribution embedded in a somewhat more extendeddark halo. We consider both spherical and non-spherical stel-lar profiles, and the dwarf’s halo was kept (nearly) sphericalin all cases. The models vary in the degree of flattening,balance between rotation and dispersion, relative contribu-tion of stars and dark matter to the total mass, and densityprofiles of the dark halo (cored or cuspy). The spherically-averaged stellar density is roughly the same in all of themand follows approximately an exponential or a King pro-file. This specific choice of the functional form has little impact on the density profile of the remnant after a cou-ple of pericentre passages, when most of the mass has beenstripped and the remaining one redistributed in response totidal torques.Spherical isotropic models are constructed using the Ed-dington inversion formula, while their flattened analoguesare created with the iterative self-consistent modelling ap-proach described in Section 5 of Vasiliev (2019a). The mod-els are defined by distribution functions (DF) in actionspace, and the total potential corresponding to the densitygenerated by the DF is computed iteratively. The propertiesof the model are thus determined by the choice of the DFfamily and its parameters. There are DF families suitable MNRAS000 Kinematic maps of the Sgr remnant. Coordinates are aligned with the apparent (non-reflex-corrected) motion of the objecton the sky (the motion is in the direction of increasing χ ); the true velocity of the Sgr centre (blue arrow in panel G) points towards theGalactic plane, which lies slightly beyond the top right corner. A grid of equatorial coordinates is shown in Panel F. Most panels alsodisplay the surface density, with contours logarithmically spaced by 0.4 dex (i.e., one magnitude).Panels A and B show the perspective-corrected mean PM µ (cid:48) χ , µ (cid:48) ξ (Equation 6); panels D and E – the PM dispersions σ χ , σ ξ . The first ofthese values is inflated due to a non-negligible thickness of the galaxy, and the difference between the two dispersions (panel G) can beused to estimate the “kinematic thickness” h kin (Equation 8). It increases from ∼ ∼ v los , GSR (Equation 7) and its dispersion. Panel J shows the photometric distance estimate, which correlates with the featuresseen in the mean PM µ (cid:48) χ parallel to the direction of motion due to perspective effects. Colour scales for v los , GSR and σ los match thoseof PM for an assumed distance D = 26 . unable to match all observational constraints (see Niederste-Ostholt et al. 2012; Gibbons et al. 2017), so we concentrateon two-component models with more centrally concentratedstellar distribution embedded in a somewhat more extendeddark halo. We consider both spherical and non-spherical stel-lar profiles, and the dwarf’s halo was kept (nearly) sphericalin all cases. The models vary in the degree of flattening,balance between rotation and dispersion, relative contribu-tion of stars and dark matter to the total mass, and densityprofiles of the dark halo (cored or cuspy). The spherically-averaged stellar density is roughly the same in all of themand follows approximately an exponential or a King pro-file. This specific choice of the functional form has little impact on the density profile of the remnant after a cou-ple of pericentre passages, when most of the mass has beenstripped and the remaining one redistributed in response totidal torques.Spherical isotropic models are constructed using the Ed-dington inversion formula, while their flattened analoguesare created with the iterative self-consistent modelling ap-proach described in Section 5 of Vasiliev (2019a). The mod-els are defined by distribution functions (DF) in actionspace, and the total potential corresponding to the densitygenerated by the DF is computed iteratively. The propertiesof the model are thus determined by the choice of the DFfamily and its parameters. There are DF families suitable MNRAS000 , 000–000 (0000) E. Vasiliev & V. Belokurov 15 10 5 0 55051015 [mas/yr] A [mas/yr] B v los , GSR [km/s] C [mas/yr] D [mas/yr] E - - los [km/s] F Kinematic thickness [kpc] G Photometric thickness [kpc] H Distance [kpc] J , data - , model a , data - , model b v los , data - v los , model c , data - , model d , data - , model e los , data - los , model f Figure 6. Kinematic maps for one of the more successful models. The first three rows display the same quantities as the actualobservations shown on Figure 5, while the last two rows are the residuals for the first two rows. MNRAS , 000–000 (0000) nternal kinematics of the Sagittarius dSph X [kpc] Z [ k p c ] v li n e s o f s i g h t X [kpc] Y X [kpc] major X [kpc] minor Figure 7. Internal kinematics of the Sgr remnant (shown is same model as in Figure 6, but there is little difference between modelsthat fit the observations similarly well). Coordinates are aligned with the orbital plane of the Sgr galaxy, so that its angular momentumis antialigned with the Y axis (which points away from the image plane), and X axis is aligned with the line of sight (the observer is at X (cid:39) − . Left panel: mean relative velocity of stars with respect to the centre of the Sgr remnant. Colour indicates the magnitude of the velocityin km s − , while streamlines show its direction. It is clear that essentially no part of the remnant is moving as a solid body; not onlyit tumbles counter-clockwise (in the same sense as the orbital angular momentum), but is also strongly sheared, as evidenced by theX-shaped flow lines. Remaining panels show the components of the velocity dispersion tensor indicated by colour. Centre-left panel displays the directionperpendicular to the orbit plane, and the other two panels – the two eigenvalues (major and minor) of the in-plane velocity dispersiontensor, with its orientation also shown by red ellipses in the centre-right panel. for disky stellar systems, with roughly exponentially declin-ing surface density profiles and nearly constant thickness, orfor more dispersion-supported oblate axisymmetric systemswith rather flexible radial density profiles, which may alsohave net rotation.Since we only follow the last stages of the Sgr disrup-tion, its density profile cannot extend much beyond the tidalradius r tidal . For an orbit with an apocentre radius around60 kpc, the mean density within r tidal is approximately10 M (cid:12) kpc − , and at the pericentre radius of ∼ 16 kpc thisincreases to 3 . × M (cid:12) kpc − . For realistic models, theinitial tidal radius is ∼ 10 kpc; it drops to ∼ ∼ − 60 km s − , and the total mass lies inthe range (1 − × M (cid:12) . The fitting strategy involves several steps. Each choice ofthe initial model still leaves room for rescaling its mass andradius. We pick up several values for the mass normaliza-tion, and for each mass determine the length scaling factor that produce a final result resembling the actual Sgr rem-nant, by running a grid of reduced-resolution simulations.At this stage, the most important comparison criteria arethe remnant mass (stellar and dark), PM and v los disper-sions, and to some extent the shape; these are all linked andless sensitive to the accuracy of the final phase-space coor-dinates. Then we iteratively adjust the initial conditions, asdescribed in the Appendix A, running full-resolution simula-tions to match the present-day position and velocity. We alsomake small adjustments to the length scale of the model, toimprove the fit for the velocity dispersions. In the end, thereis no single numerical criterion describing the fit quality, andthe process involves a lot of subjective “holistic” judgementand manual labour. In total, we considered more than a hun-dred of models, of which only a small fraction were able tosatisfy all available constraints even approximately.The PM component perpendicular to the apparent mo-tion, µ (cid:48) ξ , is insensitive to the perspective effects, and all mod-els produce very similar µ (cid:48) ξ maps, which also match observa-tions rather well. On the other hand, the parallel component µ (cid:48) χ is very sensitive to the distance gradient, hence provid-ing strong constraints on the orientation and extent of theelongated remnant. The angle between its major axis andthe orbit needs to be around 45 ◦ over a distance ∼ MNRAS000 60 km s − , and the total mass lies inthe range (1 − × M (cid:12) . The fitting strategy involves several steps. Each choice ofthe initial model still leaves room for rescaling its mass andradius. We pick up several values for the mass normaliza-tion, and for each mass determine the length scaling factor that produce a final result resembling the actual Sgr rem-nant, by running a grid of reduced-resolution simulations.At this stage, the most important comparison criteria arethe remnant mass (stellar and dark), PM and v los disper-sions, and to some extent the shape; these are all linked andless sensitive to the accuracy of the final phase-space coor-dinates. Then we iteratively adjust the initial conditions, asdescribed in the Appendix A, running full-resolution simula-tions to match the present-day position and velocity. We alsomake small adjustments to the length scale of the model, toimprove the fit for the velocity dispersions. In the end, thereis no single numerical criterion describing the fit quality, andthe process involves a lot of subjective “holistic” judgementand manual labour. In total, we considered more than a hun-dred of models, of which only a small fraction were able tosatisfy all available constraints even approximately.The PM component perpendicular to the apparent mo-tion, µ (cid:48) ξ , is insensitive to the perspective effects, and all mod-els produce very similar µ (cid:48) ξ maps, which also match observa-tions rather well. On the other hand, the parallel component µ (cid:48) χ is very sensitive to the distance gradient, hence provid-ing strong constraints on the orientation and extent of theelongated remnant. The angle between its major axis andthe orbit needs to be around 45 ◦ over a distance ∼ MNRAS000 , 000–000 (0000) E. Vasiliev & V. Belokurov Figure 6 shows the kinematic maps for one of the moresuccessful models plotted in the same way as the real obser-vations (Figure 5), and additionally their residuals (differ-ences from observed values). This model, and a number ofother similarly looking models, is able to match qualitativelythe features seen in the mean PM and v los maps, and repro-duces the dispersions reasonably well. The region of lower µ (cid:48) χ extending up to ∼ ◦ from the centre towards the trail-ing arm (Panel A) is reproduced by the model, although notacross the entire minor axis. The end of this region corre-sponds to the transition to the unbound and un-twisted partof the stream, which indicates that the 3d geometry of theremnant is represented adequately. The other PM compo-nent µ (cid:48) ξ and the line-of-sight velocity are not contaminatedby perspective effects, and the mild gradient in the residualmaps indicates some deficiencies in the fit, which may bepartly alleviated by considering kinematically more compli-cated initial conditions (this particular model was initiallyspherical and non-rotating), or by adjusting the MW poten-tial. A robust conclusion from the analysis of a large suiteof models is that the Sgr core remained a bound stellar sys-tem until very recently. Models that were too large in size ornot massive enough to withstand the tidal shocks at earlierpericenter passages produced a final configuration that wastoo elongated and more closely aligned with the orbit (i.e.without a distinct “twist” in the distance gradient). An ex-ample of such configuration is given in Figure A1; the poorfit to µ (cid:48) χ is evident. Conversely, if the initial configuration istoo tightly bound and loses only relatively little mass, thefinal state is too spherical and either has too large velocitydispersion, or is too compact to match the observed extentof the remnant along its orbit (this is again most evident inthe poor fit to µ (cid:48) χ ). For instance, this was the case for theLaw & Majewski 2010 N -body model, which we also ana-lyzed in the same way as our simulations. Figure A2 showsthe kinematic maps of a system similar to the Law & Ma-jewski 2010 model, but with slightly different orbital initialconditions tuned to match the present-day position and ve-locity of Sgr remnant; their original model has very similarfeatures.Models with initially flattened and rotating stellar dis-tribution may have several observable kinematic features.First, the gradient in the velocity components representedby µ (cid:48) ξ and v los , GSR needs not be aligned with the major axis.However, there is little evidence for such a gradient in theobservations (panels B and C in Figure 5), putting an upperlimit of a few km s − on the amount of rotation about thephotometric major axis. Second, there may be a gradientin the distance along the minor axis, which would manifestitself in the µ (cid:48) χ map; again, the data do not demonstratesuch a gradient (panel A). Third, models in which the stel-lar distribution was rapidly rotating in the same sense asthe orbital angular momentum develop a tidally inducedbar upon passing the pericentre of their orbits ((cid:32)Lokas et al.2010). These models can also match most of the featuresin the data (e.g., the orientation of the bar with respect tothe orbital velocity), but they tend to have lower dispersionsin both PM components, and moreover, they have a signif-icant residual rotation about the photometric minor axis,manifested as a non-monotonic profile of v los , GSR along themajor axis (see, e.g., figure 14 in (cid:32)Lokas et al. 2010), in dis- agreement with the observations. An example of such modelis given in Figure A3. In our experiments, models with aninitial flattening greater than 1 . ∼ 10 km s − per kpc (left panel).Even though only a fraction of this mean velocity is di-rected radially, it is still clear that the system is expandingrapidly, and thus is far from a steady state. This impliesthat any analysis method based on the equilibrium assump-tion (such as Jeans equations) would give incorrect resultsregarding the mass distribution, and N -body simulations re-main the only viable modelling approach. In addition, thestreamlines of the mean velocity are not circles or ellipses,as would be in the case of a rotating system, but rather havea characteristic X shape (contraction in one direction andexpansion in the perpendicular one). The velocity ellipsoidis very anisotropic, with the largest dispersion being in thedirection perpendicular to the orbit plane (centre-left panel)closely followed by the dispersion roughly perpendicular tothe orbital velocity (centre-right panel), and the remainingcomponent being much smaller (right panel). This again de-fies expectations for a prolate stellar system, in which thedispersion along the major axis should be larger than inother directions – and yet here it is exactly the opposite(note the orientations of the velocity ellipsoids in the thirdpanel). Unfortunately, it is nearly impossible to measure di-rectly such a small dispersion in the direction parallel to theorbital motion, because, as discussed earlier, the observedPM dispersion is dominated by the spread in distances, notin the intrinsic velocity.Figure 8 illustrates the structural properties of a bunchof models with reasonably good fit quality. The stellar andtotal mass distribution is fairly similar between all thesemodels, shown by the circular velocity curve (left panel),which measures the spherically-averaged enclosed mass pro-file. It peaks around 2 . − (cid:46) 20 km s − ,and the stellar distribution is even more concentrated (starsdominate the total density within 1 kpc in most of thesemodels). For reference, we also show the circular velocitycorresponding to the mean density of ∼ . × M (cid:12) kpc − ,which is the tidal limit at the present location of the Sgrremnant. The fact that it lies just above all the models con-firms the fact that the Sgr core is tidally disturbed down toits very centre (except the nuclear star cluster M 54, whichwe did not simulate). Centre panel shows surface densityprofiles along the major and minor axes, which match theobservations fairly well.Since the transition between the remnant and the trail-ing arm of the stream occurs around 12 − ◦ from the cen-tre, as indicated both photometrically and kinematically, wetake the enclosed mass within a fiducial radius 5 kpc as our MNRAS , 000–000 (0000) nternal kinematics of the Sagittarius dSph Radius [kpc] C i r c u l a r v e l o c i t y [ k m / s ] t o t a l s t a r s Radius [kpc] S u r f a c e d e n s i t y [ M fl / k p c ] m a j o r a x i s m i n o r a x i s Time [Gyr] M a ss w i t h i n k p c [ M fl ] totalstars Figure 8. Properties of some of the more successful models for the Sgr remnant (each model shown by a different colour). Left panel: circular velocity curve v circ ( r ) ≡ (cid:112) G M ( < r ) /r corresponding to the total enclosed mass (solid lines) and stellar mass only(dashed lines). Black dotted line shows the tidal limit at the current position, indicating that the entire Sgr remnant is tidally disturbed(equivalently, its tidal radius is zero). Centre panel: surface density profiles along the major (solid lines) and minor (dashed lines) axes. For comparison, the actual observationsare plotted in black dots. Right panel: time evolution of the total (solid) and stellar (dashed lines) mass enclosed within 5 kpc. The evolution is started some2.5 Gyr before present at an apocentre; each successive pericentre passage leads to a tidal shock and causes a sudden mass loss, andthe most recent passage actually initiates a complete unbinding of the remnant, which, nevetheless, will dissolve only gradually over thenext Gyr. mass estimate. For most successful models, the total masswithin this radius is ∼ (4 ± × M (cid:12) , of which the starscontribute around 10 M (cid:12) , in accordance with our photo-metric estimate (Section 4). The mass is mainly constrainedby the velocity dispersion, but also produces configurationsof an appropriate spatial extent and elongation. We findthat the total mass profile must be more extended than thestellar component, disfavouring mass-follows-light models.Only in this case the transition between the remnant andthe stream occurs at large enough distances without produc-ing an excessively high velocity dispersion in the centre. Forinstance, in the model shown in Figures 6, 7, stars initiallyhad a King profile with M = 1 . × M (cid:12) , R = 0 . 65 kpc,and W = 4, while the dark halo had a Gaussian profile with M = 2 . × M (cid:12) and R = 2 kpc.Right panel of Figure 8 shows the time evolution ofthe enclosed mass within a fiducial radius 5 kpc (both stel-lar and total). Models started with rather different initialmasses and concentrations all converge around the present-day mass of ∼ × M (cid:12) within this radius, of which starscontribute about a quarter. The ratio between stellar andtotal mass increases with time, and was (cid:46) 10% at the begin-ning of simulation (models with a lower initial dark matterfraction were unable to maintain a sufficiently extended rem-nant without violating the velocity dispersion constraints).What is more interesting, though, is that when we con-tinue the N -body simulation into the future, all modelsdemonstrate very similar behaviour: the mass within a fixedradius drops precipitously – in other words, the Sgr galaxyis completely disrupted over its next orbit. Although someconcentration of mass close to its centre-of-mass is retained,it no longer remains a gravitationally bound system. Thuswe conclude that we are witnessing the final demise of thisonce third-massive satellite of our Galaxy. Despite its proximity, the dynamical state of the Sgr dSphhas been studied only in a few papers. The most widelyknown N -body model of Law & Majewski (2010) focused onreproducing the properties of the Sgr stream, not its core. Asmentioned earlier, the remnant is much too compact in thismodel, and even though its line-of-sight velocity dispersiondoes not exceed the observational limits, the PM dispersionis too low, and the extent of the region of low µ (cid:48) χ seen inthe data is not reproduced by the model, since the transitionfrom the nearly spherical core to the stream occurs too early(Panels A, G, J in Figure A2). They estimate the total massof the remnant to be 2 . +1 . − . × M (cid:12) based on the velocitydispersion in the stream; the actual remnant mass in the N -body model is not quoted, but is likely smaller than ourestimate ( ∼ × M (cid:12) ) based on the lower velocity dis-persion. The smooth monotonic trend of v los , GSR along themajor axis, seen in the data, is not quite reproduced by themodel: it displays a “kink” (sudden steepening of the gradi-ent) in the v los , GSR within 2 − ◦ from the centre, where theremnant is intrinsically not rotating (Panel C in the abovefigure). Based on the analysis of our suite of models, thisfeature is characteristic of more concentrated systems thatresist the tidal perturbation in their central parts.Pe˜narrubia et al. (2010) proposed a model in which theSgr galaxy was initially a rapidly rotating disc embedded in adark halo. This scenario could explain the bifurcation in theSgr stream, but predicted a high degree of rotation in theremnant, which subsequently was not confirmed by newerobservations in Pe˜narrubia et al. (2011). Notably, the modelpredicted a strong gradient in v los , GSR along both axes (ma-jor and minor), whereas in reality the variation along themajor axis has the opposite sign and is much shallower, asevidenced by larger-scale kinematic measurements of Frinch-aboy et al. (2012); the minor axis gradient in the data has MNRAS000 MNRAS000 , 000–000 (0000) E. Vasiliev & V. Belokurov the same sign but is also much weaker than in the model.We also analyzed the PM field predicted by this model andfound it to be similarly discrepant with the observationaldata.(cid:32)Lokas et al. (2010) presented another scenario with adisky Sgr progenitor, whose internal angular momentum wasinitially nearly co-aligned with its orbital angular momen-tum (i.e., the rotation is prograde with respect to the orbit).In this case, a strong bar perturbation is induced during apericentre passage, and they find that after a second pas-sage, their model provides a good match to the observations(in particular, the mean v los and its dispersion measured byFrinchaboy et al. 2012). A more detailed analysis of theirfigures 10 and 11 suggests that the remnant is still too hotkinematically ( σ los ∼ − 20 km s − as opposed to the ob-served values 12 − 15 km s − ), and the v los , GSR profile alongthe major axis (figure 14) is non-monotonic, unlike the ob-served one, although the disagreement is at the level of only afew km s − . The enclosed mass within 5 kpc is ∼ × M (cid:12) ,similar to our estimates, but it is more centrally concen-trated (the peak circular velocity of ∼ 21 km s − is attainedat ∼ . . − v los , GSR along the major axis. While we do not rule out apossibility of an initially rotating Sgr progenitor, it appearsto be unnecessary to explain the shape and kinematics of itsremnant (although still might be needed to reproduce somefeatures of the Sgr stream). The N -body models considered in this paper were specifi-cally designed to reproduce the properties of the Sgr rem-nant, not the stream, unlike most existing studies. Our treat-ment of the Sgr galaxy orbit is rather simplistic: all thatwe require is to arrive at the present location with the ve-locity consistent with observations. We do follow the Sgrgalaxy as a live N -body system, however, the MW is rep-resented as a static external potential, we ignore the effectsof the LMC, and limit ourselves to one particular choice ofthe MW potential. We also neglect the dynamical friction,which is likely unimportant over the last 1 Gyr owing to arelatively small mass of the remnant, but certainly plays arole at earlier times.We also start our simulations relatively recently (2 − . . − 2) after the first pericentre passage, but much more (factor of 2 − 5) after the second (penultimate)one. Our models do not necessarily represent the mass losshistory particularly well at early times, but they should bemore reliable over the last orbital period and can be usedto forecast the future evolution, being well constrained bythe present-day state of the Sgr remnant. In short, the 2Gyr of evolution is just a device to produce a realisticallylooking tidally perturbed model, which is then comparedwith observations. Our approach is also driven by practicalconsiderations: we find it vitally important to obtain a veryaccurate fit to the present-day position and velocity in orderto adequately compare the kinematics of different models.A longer simulation period and incorporation of additionalperturbations on the orbit will make this task still harder.The question remains whether the parameters responsi-ble for the details of the preceding evolution (i.e., potential,orbit shape, mass loss history) can leave a discernible sig-nature in the present-day state of the remnant. Based onour preliminary experiments, the answer is likely positive,and the models can and should be refined while adapting toother observational constraints on the structure of the Sgrstream. However, we believe that the main features of theSgr remnant are unlikely to change qualitatively. We presented a detailed analysis of the Sgr galaxy core, us-ing the data from the Gaia DR2 catalogue and other ex-isting spectroscopic datasets. We developed a multidimen-sional mixture model to classify the catalogue into candidateSgr members and field stars, which simultaneously producesastrometric kinematic maps of the Sgr galaxy. The finallist of candidate members contains ∼ . × stars withextinction-corrected G -band magnitude brighter than 18 (upto and including the red clump). The observational datasummarized on Figure 5 already reveal a number of impor-tant features pertaining to both the spatial and kinematicproperties of the Sgr remnant and its transition into the Sgrstream. A more detailed interpretation was conducted withthe help of a large suite of N -body models for the Sgr rem-nant, which were evolved in the Galactic tidal field for 2.5orbital periods preceding the most recent pericentre passage(which occurred ∼ 30 Myr ago). The results of this analysiscan be summarized as follows. • We estimated the total stellar mass of the Sgr remnantto be ∼ M (cid:12) from its photometry. • At the same time, the total mass of the Sgr remnantwithin 5 kpc from its centre is a factor of four higher. Wefind that a mass-follows-light model is a worse match for theobservational constraints compared to models with a moreextended dark halo than the stellar distribution. • The 3d shape and orientation of the Sgr remnant isstrongly constrained by the imprint it leaves in the PM field.We find that the remnant is a prolate structure tilted at ∼ ◦ with respect to its orbit, which transitions into thetidal stream beyond ∼ • The observed cigar-like shape is caused by the Galactictidal field and is well reproduced even by models that start MNRAS , 000–000 (0000) nternal kinematics of the Sagittarius dSph as spherically symmetric, without the need to invoke initialflattening or rotation. • The combination of a relatively low velocity dispersionwith an extended prolate shape strongly suggests that theSgr remnant ceased to be gravitationally bound after themost recent pericentre passage, and will gradually dissolveover the next orbital period. • The remnant is significantly out of equilibrium to renderthe classical dynamical modelling methods useless.Our models were specifically tailored to reproduce theobserved properties of the Sgr remnant, not the stream, andare a poor fit to the latter. It is likely that a successful recon-struction of the stream would require to add more physicalingredients influencing the Sgr orbit, which could also af-fect the properties of the remnant. Nevertheless, its globalfeatures are tightly constrained by abundant observationaldata, and should be taken into account by any study focus-ing on the Sgr stream or on the perturbations in the MWdisc produced by the Sgr galaxy.We thank Jorge Pe˜narrubia and Denis Erkal for pro-viding snapshots of their simulations and for enlighten-ing discussions. We used the Whole Sky Database (wsdb)created and maintained by Sergei Koposov with financialsupport from the Science & Technology Facilities Coun-cil (STFC). This work uses the data from the EuropeanSpace Agency mission Gaia ( ), processed by the Gaia Data Processing andAnalysis Consortium ( ). This work started at the Kavli In-stitute of Theoretical Physics during the programme “Dy-namical Models for Stars and Gas in Galaxies in the GaiaEra”, which was supported in part by the National ScienceFoundation under Grant No. NSF PHY-1748958. REFERENCES Alcock C., Allsman R., Alves D., et al., 1997, ApJ, 474, 217Alfaro-Cuello M., Kacharov V., Neumayer N., et al., 2020, ApJ,892, 20Antoja T., Helmi A., Romero-G´omez M., et al., 2018, Natur, 561,360Antoja T., Ramos P., Mateu C., Helmi A., Anders F., Jordi C.,Carballo-Bello J., 2020, A&A, 635, L3Astropy collaboration (Price-Whelan et al.), 2018, ApJ, 156, 123Bailin J., 2003, ApJL, 583, L79Baumgardt H., Hilker M., Sollima A., Bellini A., 2019, MNRAS,482, 5138Bellazzini M., Ibata R., Chapman S., et al., 2008, AJ, 136, 1147Belokurov V., Zucker D., Evans N.W., et al., 2006, ApJ, 642, L137Bland-Hawthorn J., Sharma S., Tepper-Garcia T., et al., 2019,MNRAS, 486, 1167Bovy J., Hogg D., Roweis S., 2011, Ann. Appl. Stat., 5, 1657Cappellari M., Copin Y., 2003, MNRAS, 342, 345Carrillo I., Minchev I., Steinmetz M., et al., 2019, MNRAS, 490,797Cseresnjes P., Alard C., Guibert J., 2000, A&A, 357, 871Cseresnjes P., 2001, A&A, 375, 909Dehnen W., 2000, ApJ, 536, L9; ascl:1402.031Dempster A., Laird N., Rubin D., 1977, J. R. Stat. Soc. series BMethodological, 39, 1Erkal D., Belokurov V., Laporte C., et al., 2019, MNRAS, 487,2685 Ferguson P., Strigari L., 2020, MNRAS (in press);arXiv:1909.11103Frinchaboy P., Majewski S., Mu˜noz R., et al., 2012, ApJ, 756, 74 Gaia Collaboration (Brown et al.), 2018a, A&A, 616, 1 Gaia Collaboration (Babusiaux et al.), 2018b, A&A, 616, 10Garavito-Camargo N., Besla G., Laporte C., Johnston K. , G´omezF., Watkins L., 2019, ApJ, 884, 51Gibbons S., Belokurov V., Evans N. W., 2017, MNRAS, 464, 794G´omez F., Besla G., Carpintero D., Villalobos ´A., O’Shea B., BellE., 2015, ApJ, 802, 128Hamanowicz A., Pietrukowicz P., Udalski A., et al., 2016, ActaAstronomica, 66, 197Helmi A., 2004, ApJL, 610, L97Ibata R., Gilmore G., Irwin M., 1994, Nature, 370, 194Ibata R., Gilmore G., Irwin M., 1995, MNRAS, 277, 781Ibata R., Wyse R., Gilmore G., Irwin M., Suntzeff N., 1997, AJ,113, 634Ibata R., Razoumov A., 1998, A&A, 336, 130Ibata R., Bellazzini M., Thomas G., Malhan K., Martin N.,Famaey B., Siebert A., 2020, ApJ, 891, L19Jiang I.-G., Binney J., 2000, MNRAS, 314, 468Johnston K., Law D., Majewski S., 2005, ApJ, 619, 800Kunder A., Chaboyer B., 2009, AJ, 137, 4478Laporte C., Johnston K., G´omez F., Garavito-Camargo N., BeslaG., 2018, MNRAS, 481, 286Laporte C., Minchev I., Johnston K., G´omez F., 2019, MNRAS,485, 3134Law D., Majewski S., 2010, ApJ, 714, 229Lindegren L., Hernandez J., Bombrun A., et al., 2018, A&A, 616,2(cid:32)Lokas E., Kazantzidis S., Majewski S., Law D., Mayer L., Frinch-aboy P., 2010, ApJ, 725, 1516Majewski S., Skrutskie M., Weinberg M., Ostheimer J., 2003,ApJ, 599, 1082Majewski S., Hasselquist S., (cid:32)Lokas E., et al., 2013, ApJ, 777, L13Marconi G., Buonanno R., Castellani M., Iannicola G., MolaroP., Pasquini L., Pulone L., 1998, A&A, 330, 453Mateo M., Udalski A., Szymanski M., Kaluzny J., Kubiak M.,Krzeminski W., 1995, AJ, 109, 588Mateo M., Olszewski E. W., Morrison H. L., 1998, ApJL, 508,L55McMillan P., 2017, MNRAS, 465, 76Niederste-Ostholt M., Belokurov V., Evans N. W., Pe˜narrubia J.,2010, ApJ, 712, 516Niederste-Ostholt M., Belokurov V., Evans N. W., 2012, MNRAS,422, 207Pe˜narrubia J., Belokurov V., Evans N. W., et al., 2010, MNRAS,408, L26Pe˜narrubia J., Zucker D., Irwin M., et al., 2011, ApJ, 727, L2Purcell C., Bullock J., Tollerud E., Rocha M., Chakrabarti S.,2011, Nature, 477, 301Ramos P., Mateu C., Antoja T., Helmi A., Castro-Ginard A.,Balbinot E., 2020, arXiv:2002.11142Teuben P., 1995, in Shaw R. A., Payne H. E., Hayes J. J. E., eds,ASP Conf. Ser. Vol. 77, Astronomical Data Analysis Softwareand Systems IV. Astron. Soc. Pac., San Francisco, p. 398;ascl.net/1010.051Vasiliev E., 2019a, MNRAS, 482, 1525Vasiliev E., 2019b, MNRAS, 489, 623Vera-Ciro C., Helmi A., 2013, ApJ, 773, L4 APPENDIX A: DETERMINATION OF THEORBIT INITIAL CONDITIONS The N -body simulations discussed in Section 6 start ∼ . MNRAS , 000–000 (0000) E. Vasiliev & V. Belokurov present-day centre-of-mass position and velocity of the Sgrremnant, using the procedure detailed below.The first guess for the initial conditions (IC) of theSgr orbit comes from integrating the orbit backwards in thestatic MW potential; however, the actual trajectory of a dis-rupting satellite deviates from a test-particle orbit, necessi-tating further refinement of the initial conditions. We em-ploy the standard Gauss–Newton iterative procedure to findthe orbital IC w in ≡ { x , v } in leading to the given final po-sition and velocity w end after a fixed time T end . For a givenchoice of w in , , we follow an ensemble of simulations withslightly offset IC w in ,k , k = 1 ..K , which produce the endstates w end ,k . Then the Jacobian matrix J ≡ ∂ w end /∂ w in isapproximated by finite differences: J ≈ δ w end δ w − , wherethe columns of matrices δ w ... contain the difference vectors w ...,k − w ..., . Finally, the next choice of IC is given by w (new)in , = w in , − J − ( w end , − w trueend ).The procedure outlined above is a general way of solv-ing nonlinear equation systems, but to make it practical inthe present case, a few adaptations were made. Since theSgr remnant has just passed the pericentre, its position andvelocity are rapidly varying. The main effect of perturbingthe IC is the slight change in the orbital energy and thecorresponding change in flight time. After 2 Gyr of evolu-tion, this translates to the final states of perturbed orbitsbeing stretched along the trajectory, making the Jacobianextremely degenerate, with condition number exceeding 10 .On the other hand, this orbital motion is fairly predictableand may be treated separately from the perturbations per-pendicular to the orbit.We introduce two auxiliary coordinate systems for theinitial and final states, aligned with the position and ve-locity vectors of the unperturbed orbit. Namely, w ( t ) = w in , + B in p ( t ) = w trueend + B end q ( t ), where p ( t ) and q ( t )are the phase-space coordinates of an orbit in either ofthe two auxiliary systems, and the orthogonal matrices B in , B end are defined in such a way that the orbital mo-tion occurs along the first component of p or q at t = 0or t = T end , respectively. The first column of B in is theunit-normalized time derivative of the 6d phase-space co-ordinates, i.e., { v in , , − ∂ Φ /∂ x | x = x in , } , and the remainingcolumns are all orthogonal to the first column and betweenthemselves, but otherwise arbitrary. Similarly, B end is de-fined by the velocity and acceleration at w trueend .The IC of the baseline orbit is thus p in , ≡ p ( t =0) = , and the ICs of K perturbed orbits p in ,k , k = 1 ..K are confined to the 5d subspace defined by setting the firstcomponent p (1) k of each vector to zero, i.e., orthogonal tothe unperturbed orbit at the initial point. The final states q end ,k ≡ q k ( T end ) of all K + 1 orbits (including the un-perturbed one) do not necessarily have zero in their firstcomponent (of course, the goal is to have q end , = , but weare searching for this solution iteratively). Nevertheless, allorbits in the bundle do cross the reference subspace q (1) = 0at some moment of time (which may be greater or less than Since the position and velocity have different units, we intro-duce dimensional scaling factors before rotating the 6d phasespace, choosing them to approximately match the magnitude ofposition and velocity variations – in this case, using 1 kpc and 10km s − as scale factors. T end , hence we run the simulation for a slightly longer timeto ensure the crossing of this subspace). We determine thetime of crossing and the corresponding coordinates for each k -th orbit in the following way. In the N -body simulation, westore the position and velocity of the Sgr remnant’s centre ateach timestep w k ( t ) and linearly transform it to q k ( t ). Thesevalues are still subject to numerical noise, so we first locatethe snapshot closest to q (1) k ( t ) = 0, and then fit a smooth or-bit to the centre positions over the interval ± . T cross ,k and the corresponding 5d coordinates q (2 − ,k .The key point now is that instead of using the lin-early transformed initial and end states p in ,k , q end ,k inthe Jacobian (which, of course, would not change its con-dition number), we define the end state ˜ q end ,k of each orbitby the time of crossing and the corresponding 5d coordi-nates. Of course, the correspondence between q end ,k and˜ q end ,k ≡ { T cross ,k , q (2 − ,k } is well defined, however, thistransformation is nonlinear in a favourable way. The mainnonlinear effect is the motion along a curved trajectory,which can be followed explicitly, e.g., by numerically inte-grating a test-particle trajectory in the static MW potential(over such short timescales the actual trajectory of the Sgrremnant’s centre in the simulation is well approximated by atest-particle orbit). Similarly, the initial states of actual or-bits are already confined to a 5d subspace, but for any pointoutside this subspace, it can be projected back by followinga curvilinear trajectory until crossing p (1) = 0, thus defin-ing the nonlinearly transformed vector ˜ p that contains the“flight time” in the first component and the 5d coordinatesof the crossing point in the remaining components. Oncethe curved orbital motion is compensated, the Jacobian oftransformation between the initial ˜ p in ,k and the final ˜ q end ,k states is much better behaved, with the condition number (cid:46) . 05 kpc and 0 . − in 3 − MNRAS , 000–000 (0000) nternal kinematics of the Sagittarius dSph 15 10 5 0 55051015 [mas/yr] A [mas/yr] B v los , GSR [km/s] C [mas/yr] D [mas/yr] E - - los [km/s] F Z v [km/s] li n e s o f s i g h t G Photometric thickness [kpc] H Distance [kpc] J , data - , model a , data - , model b v los , data - v los , model c , data - , model d , data - , model e los , data - los , model f Figure A1. Same as Figure 6, but for a model that is too much tidally stretched and is more aligned with the stream (Panel G showingthe side-on view, as in Figure 7). The trailing side is too close (Panel J), causing a serious misfit in µ (cid:48) χ (Panel A).MNRAS000 Same as Figure 6, but for a model that is too much tidally stretched and is more aligned with the stream (Panel G showingthe side-on view, as in Figure 7). The trailing side is too close (Panel J), causing a serious misfit in µ (cid:48) χ (Panel A).MNRAS000 , 000–000 (0000) E. Vasiliev & V. Belokurov 15 10 5 0 55051015 [mas/yr] A [mas/yr] B v los , GSR [km/s] C [mas/yr] D [mas/yr] E - - los [km/s] F Z v [km/s] li n e s o f s i g h t G Photometric thickness [kpc] H Distance [kpc] J , data - , model a , data - , model b v los , data - v los , model c , data - , model d , data - , model e los , data - los , model f Figure A2. Same as Figure 6, but for the Law & Majewski (2010) model (resimulated with a slightly different orbit to better matchthe present-day position and velocity of Sgr), which is too compact and transitions to the stream too early (Panels A, G).MNRAS , 000–000 (0000) nternal kinematics of the Sagittarius dSph 15 10 5 0 55051015 [mas/yr] A [mas/yr] B v los , GSR [km/s] C [mas/yr] D [mas/yr] E - - los [km/s] F Z v [km/s] li n e s o f s i g h t G Photometric thickness [kpc] H Distance [kpc] J , data - , model a , data - , model b v los , data - v los , model c , data - , model d , data - , model e los , data - los , model f Figure A3. Same as Figure 6, but a model that is initially rotating and forms a tidally induced bar, as in the (cid:32)Lokas et al. (2010)scenario. The residual rotation is manifested in the non-monotonic v los , GSR profile along the major axis (Panel C).MNRAS000 Same as Figure 6, but a model that is initially rotating and forms a tidally induced bar, as in the (cid:32)Lokas et al. (2010)scenario. The residual rotation is manifested in the non-monotonic v los , GSR profile along the major axis (Panel C).MNRAS000 , 000–000 (0000) E. Vasiliev & V. Belokurov Table A1. Measurements of the mean PM and its dispersion in200 Voronoi bins, each bin containing ∼ stars. α and δ arethe average coordinates of stars in each bin; µ α and µ α are themean PM components, σ α and σ δ are the dispersions, and ρ isthe correlation coefficient, representing the non-diagonal elementof the 2d dispersion tensor. Statistical uncertainty on PM dis-persion is (cid:46) . 01 mas yr − , and systematic error in mean PM is (cid:46) . 05 mas yr − . α δ µ α µ δ σ α σ δ ρ deg mas yr − mas yr − Table A1 – continued Mean PM and its dispersion α δ µ α µ δ σ α σ δ ρ , 000–000 (0000) nternal kinematics of the Sagittarius dSph Table A1 – continued Mean PM and its dispersion α δ µ α µ δ σ α σ δ ρ Table A1 – continued Mean PM and its dispersion α δ µ α µ δ σ α σ δ ρ Table A2. Measurements of the line-of-sight velocity and its dis-persion in 36 Voronoi bins. α and δ are the average coordinatesof stars in each bin; v los is the mean heliocentric line-of-sightvelocity, σ los is its dispersion, and (cid:15) v , (cid:15) σ are their statistical un-certainties. The final column is the number of stars in the bin. α δ v los σ los (cid:15) v (cid:15) σ N (cid:63) deg km s − km s −000