The assembly of our Galaxy can be reconstructed using the motions and chemistry of individual starsSUP12SUP. Chemodynamical studies of the stellar halo near the Sun have indicated the presence of multiple componentsSUP3SUP such as streamsSUP4SUP and clumpsSUP5SUP as well as correlations between the stars chemical abundances and orbital parametersSUP68SUP. Recently analyses of two large stellar surveysSUP910SUP revealed the presence of a well populated elemental abundance sequenceSUP711SUP two distinct sequences in the colourmagnitude diagramSUP12SUP and a prominent slightly retrograde kinematic structureSUP1314SUP in the halo near the Sun which may trace an important accretion event experienced by the GalaxySUP15SUP. However the link between these observations and their implications for Galactic history is not well understood. Here we report an analysis of the kinematics chemistry age and spatial distribution of stars that are mainly linked to two major Galactic components the thick disk and the stellar halo. We demonstrate that the inner halo is dominated by debris from an object that at infall was slightly more massive than the Small Magellanic Cloud and which we refer to as GaiaEnceladus. The stars that originate in GaiaEnceladus cover nearly the full sky and their motions reveal the presence of streams and slightly retrograde and elongated trajectories. With an estimated mass ratio of four to one the merger of the Milky Way with GaiaEnceladus must have led to the dynamical heating of the precursor of the Galactic thick disk thus contributing to the formation of this component approximately ten billion years ago. These findings are in line with the results of galaxy formation simulations which predict that the inner stellar halo should be dominated by debris from only a few massive progenitorsSUP216SUP.
['Astrophysics - Astrophysics of Galaxies']
The merger that led to the formation of the Milky Ways inner stellar halo and thick disk
{"The merger that led to the formation of the Milky Way's inner stellar halo and thick disk": "Amina Helmi 1 , Carine Babusiaux 2 , Helmer H. Koppelman 1 , Davide Massari 1 , Jovan Veljanoski 1 , Anthony G. A. Brown 3  \nThe assembly process of our Galaxy can be retrieved using the motions and chemistry of individual stars. 1, 2 Chemo-dynamical studies of the nearby halo have long hinted at the presence of multiple components such as streams, 3 clumps, 4 duality 5 and correlations between the stars' chemical abundances and orbital parameters. 6, 7, 8 More recently, the analysis of two large stellar surveys 9, 10 have revealed the presence of a well-populated chemical elemental abundance sequence, 7, 11 of two distinct sequences in the colour-magnitude diagram, 12 and of a prominent slightly retrograde kinematic structure 13,14 all in the nearby halo, which may trace an important accretion event experienced by the Galaxy. 15 Here report an analysis of the kinematics, chemistry, age and spatial distribution of stars in a relatively large volume around the Sun that are mainly linked to two major Galactic components, the thick disk and the stellar halo. We demonstrate that the inner halo is dominated by debris from an object which at infall was slightly more massive than the Small Magellanic Cloud, and which we refer to as Gaia-Enceladus. The stars originating in Gaia-Enceladus cover nearly the full sky, their motions reveal the presence of streams and slightly retrograde and elongated trajectories. Hundreds of RR Lyrae stars and thirteen globular clusters following a consistent age-metallicity relation can be associated to Gaia-Enceladus on the basis of their orbits. With an estimated 4:1 mass-ratio, the merger with Gaia-Enceladus must have led to the dynamical heating of the precursor of the Galactic thick disk and therefore contributed to the formation of this component approximately 10 Gyr ago. These findings are in line with simulations of galaxy formation, which predict that the inner stellar halo should be dominated by debris from just a few massive progenitors. 2, 16  \nThe sharp view provided by the second data release (DR2) of the Gaia mission, 17 has recently revealed 14 that, besides a few tight streams, a significant fraction of the halo stars  \nnear the Sun are associated with a single large kinematic structure that has slightly retrograde mean motion and which dominates the Hertzsprung-Russell diagram's (HRD) blue sequence revealed in the Gaia data. 12 This large structure is readily apparent (in blue) in Fig. 1a, which shows the velocity distribution of stars (presumably belonging to the halo) in the Solar vicinity inside a volume of 2.5 kpc radius from Gaia data (see Methods for details). Fig. 1b shows the velocity distribution from a simulation of the formation of a thick disk via a 20% mass-ratio merger. 18 The similarity between the panels suggests that the retrograde structure could be largely made up of stars originating in an external galaxy that merged with the Milky Way in the past.  \n<!-- image -->  \nFig. 1.Velocity distribution of stars in the Solar vicinity in comparison to a merger simulation . In the left panel, the velocities of stars in the disk are plotted with grey density contours (because of the large number of stars), while the halo stars (selected as those with | v -v LSR | > 210 km/s, where v LSR is the velocity of the Local Standard of Rest) are shown as points. The blue points are part of a prominent structure with slightly retrograde mean rotational motion, and have been selected here as those having -1500 < L z < 150 kpc km/s and energy E > -1 . 8 × 10 5 km 2 /s 2 (see Methods for details). The panel on the right shows the distribution of star particles in a small volume extracted from a simulation 18 of the formation of a thick disk via a 5:1 merger between a satellite (in blue) and a pre-existing disk (in black). The overall morphology and the presence of an arch (from V y ∼ -450 km/s and V ⊥ = √ V 2 x + V 2 z ∼ 50 km/s to V y ∼ -150 km/s and V ⊥ ∼ 300 km/s seen in the left panel) can be reproduced qualitatively after appropriately scaling the velocities (see Methods), in a simulation where the satellite is disky (rather than spherical, as the arch-like feature is sharper), and on a retrograde orbit inclined by ∼ 30 o to 60 o .  \n<!-- image -->  \nSupport for this hypothesis comes from the chemical abundances of stars provided by the APOGEE survey. 9 In Fig. 2a we plot the [ α /Fe] vs [Fe/H] abundances for a sample of stars cross-matched to Gaia DR2 (see Methods for details). α -elements are produced by massive  \nstars that die fast as supernovae (SNII), while iron, Fe, is also produced in SNI explosions of binary stars. Therefore in a galaxy, [ α /Fe] decreases with time (as [Fe/H] increases). Fig. 2a shows the well-known sequences defined by the thin and thick disks. The vast majority of the retrograde structure's stars (in blue), follow a well-defined separate sequence that extends from low to relatively high [Fe/H]. (The presence of lowα stars with retrograde motions in the nearby halo has in fact been reported before 7, 19 but for a small sample. The existence of a well-populated sequence with lower [ α /Fe] was demonstrated very recently using also APOGEE data 11 ). An independent analysis 15 has confirmed the relation between Gaia 's HRD blue sequence and the kinematic structure shown in Fig. 1a, and established firmly the link to the lowα stars using both earlier data 7 as well as APOGEE, thereby putting the accretion hypothesis on more secure ground.  \nThe large metallicity spread of the retrograde structure stars depicted in Fig. 2b, implies that they did not form in a single burst in a low mass system. Furthermore because the more metal-rich stars have lower [ α /Fe] at the characteristic metallicity of the thick disk ([Fe/H] ∼ -0 . 6), this means that they were born in a system with a lower star formation rate than the thick disk. The star formation rate required to match the α -poor sequence of the APOGEE data has recently been calculated using a chemical evolution model and including different elemental abundances, 20 and found to be 0 . 3 M glyph[circledot] /yr lasting for about 2 Gyr. This implies a stellar mass for the progenitor system of ∼ 6 × 10 8 M glyph[circledot] , a value that is consistent with the large fraction of nearby halo stars being associated with the structure given estimates of the local halo density, 21 and which is comparable to the present-day mass of the Small Magellanic Cloud. 22 Interestingly, previous work 11 has shown that the trends in the abundances of low metallicity stars in the Large Magellanic Cloud actually overlap quite well with the sequence, implying that the structure was comparable to the Large Magellanic Cloud in its early years. Furthermore and perhaps even more importantly, because [ α /Fe] must decrease as [Fe/H] increases, the stars in the structure could not have formed in the same system as the vast majority of stars in the Galactic thick disk. They must have formed, as previously suggested, 7, 14, 15 in a separate galaxy, which we refer to as Gaia-Enceladus hereafter (see Methods for the motivation behind the naming).  \nWe now explore whether the Gaia-Enceladus galaxy could have been responsible at least partly for the formation of the thick disk, 13,14,15 as the comparison between the data and the simulation shown in Fig. 1 would suggest. In that case, a pre-existing disk must have been in place at the time of the merger. Fig. 2c plots the HRD of the halo stars in Fig. 1a showing the Gaia-Enceladus stars (in blue) populating Gaia 's blue sequence. 12,14,15 The thinness of this sequence is compatible with an age range from ∼ 10 to 13 Gyr given the stars' abundance sequence, as indicated by the plotted isochrones. 23 Previous studies, 24,25 on which this age range is based, have shown that the stars on the α -poor sequence are  \nFig. 2.Astrophysical properties of stars in Gaia-Enceladus. Panel a) shows the chemical abundances for a sample of stars located within 5 kpc from the Sun resulting from the cross-match between Gaia and APOGEE. The blue circles correspond to 590 stars that have -1500 < L z < 150 km/s kpc and E > -1 . 8 × 10 5 km 2 /s 2 (as in Fig. 1a, but now for a larger volume to increase the sample size, see Methods). Note the clear separation between the thick disk and the sequence defined by the majority of the stars in the retrograde structure, except for a small amount of contamination (17%) by thick disk stars (i.e. on the α -rich sequence) that share a similar phasespace distribution as the structure. The error bar in the lower left corner shows the median error for the sample. The solid (dotted) histogram in panel b) shows the metallicity distribution of the structure without (with) the subset of α -rich stars. Their distribution peaking at [Fe/H] ∼ -1 . 6, is very reminiscent of that of the stellar halo. 21 Panel c) is the HRD for halo stars (black points, selected as in Fig. 1a with the additional photometric quality cuts: 12 E ( B -V ) < 0 . 015 to limit the impact in the magnitudes and colours to less than 0.05 mag, and phot-bp-rp-excess-factor < 1 . 3 + 0 . 06( G BP -G RP ) 2 ) and shows Gaia 's blue and red sequences. Gaia-Enceladus stars are plotted with dark blue symbols, with those in APOGEE within 5 kpc and with [ α /Fe] < -0 . 14 -0 . 35 [Fe/H], in light blue. The superimposed isochrones 23 based on previous work 25 show that an age range from 10 to 13 Gyr is compatible with the HRD of Gaia-Enceladus.  \n<!-- image -->  \n<!-- image -->  \nM  \nyounger than those on the α -rich sequence for -1 < [Fe/H] < -0 . 5. This implies that the progenitor of the Galactic thick disk was in place when Gaia-Enceladus fell in, which based on the ages of its youngest stars, would suggest that the merger took place around 10 Gyr ago, i.e. at redshift z ∼ 1 . 8.  \nSuch a prominent merger must have left debris over a large volume of the Galaxy. To explore where we may find other tentative members of Gaia-Enceladus beyond the solar neighbourhood, we consider stars in the Gaia 6D sample with 20% relative parallax error, with glyph[pi1] > 0 . 1 mas and having -1500 < L z < 150 kpc km/s. Fig. 3 shows that nearby tentative Gaia-Enceladus stars (with glyph[pi1] > 0 . 25 mas, darker points) are distributed over the  \nFig. 3.Sky distribution of tentative Gaia-Enceladus members from a Gaia subsample of stars with full phase-space information . These stars have glyph[pi1] > 0 . 1 mas, relative parallax error of 20%, and are colour-coded by their distance from the Sun (from near in dark red to far in light yellow). They satisfy the condition -1500 < L z < 150 kpc km/s. Because of the larger volume explored, we do not include additional selection criteria based on energy, as done for Fig. 2 (since energy depends on the Galactic potential whose spatial variation across the volume explored is less well-constrained than its local value), nor on velocity as for Fig. 1a (because of spatial gradients). We thus expect some amount of contamination by thick disk stars, especially towards the inner Galaxy (see Methods). The starry symbols are Gaia RR Lyrae stars potentially associated to this structure. To identify these, we bin the sky in 128 × 128 elements, and log glyph[pi1] in bins of 0.2 width (mimicking the relative parallax error), and measure the average proper motion of Gaia-Enceladus stars in each 3D bin. We then require that the RR Lyrae have the same proper motion (within 25 km/s in each component at their distance), which for example corresponds to 1 mas/yr for those with glyph[pi1] ∼ 0 . 2 mas. Globular clusters with L z < 250 kpc km/s, located between 5 and 15 kpc from the Sun, and 40 o away from the Galactic centre, are indicated with solid circles.  \n<!-- image -->  \nwhole sky, this subset being more than 90% complete. More distant stars are preferentially found in specific regions of the sky, and although for such small Gaia parallaxes ( glyph[pi1] = 0 . 1 -0 . 25 mas) the zero-point offset ( ∼ -0 . 03 mas) is significant and this affects the selection in L z , it does not to the extent that it can produce the observed asymmetry on the sky. At least in part this asymmetry is due to the 20% relative parallax error cut, as highlighted in Fig. 4 (see Methods for more details and also for possible links to known overdensities). In Figure 3 we have also overplotted (with starry symbols) a subset of 200 Gaia RR Lyrae stars. 26 These have proper motions similar to the mean of the candidate Gaia-Enceladus stars with full phase-space information, at their sky position and parallax. Thirteen globular  \nclusters can also be associated to Gaia-Enceladus on the basis of their angular momenta 27 (Fig. 3). All these clusters show a consistent age-metallicity relation. 28  \nFig. 4.Kinematic properties of Gaia-Enceladus tentative members on the sky. The plotted stars are a subset of those in Fig. 3 with 0 . 1 < glyph[pi1] < 0 . 2 mas, and are colour coded by their radial velocity with the arrows indicating their direction of motion. To avoid cluttering the panels correspond to different proper motion ranges, and we have removed stars close to the bulge (within 30 o in longitude and 20 o in latitude). All velocities have been corrected for the Solar and for the Local Standard of Rest motions. The grey contours encompass 90% of the stars in the 6D Gaia set with 0 . 1 < glyph[pi1] < 0 . 2 mas and having 20% relative parallax error, and clearly demonstrate this selection criterion impacts our ability to identify distant Gaia-Enceladus stars in certain regions of the sky. Notice the large-scale pattern in the radial velocity, as well as its correlation with the proper motion component µ l : stars with µ l > 0 (top panels) have v GSR > 0 for l glyph[greaterorsimilar] 75 o and v GSR < 0 for l glyph[lessorsimilar] -75 o , while the opposite occurs for µ l < 0. Such a global pattern (and its reversal for -75 o glyph[lessorsimilar] l glyph[lessorsimilar] 75 o ), arises because of the coherent retrograde sense of rotation of the stars in their orbits (i.e. they have mostly L z glyph[lessorsimilar] 0), but the correlation with µ l is a result of their elongated orbits, e.g. we see that if µ l > 0 and l glyph[greaterorsimilar] 75 o stars are typically moving outwards with high speed and away from the Solar radius ( v GSR glyph[greaterorsimilar] 100 km/s).  \n<!-- image -->  \nFig. 4 shows the velocity field of the more distant stars associated to Gaia-Enceladus. Notice the large-scale gradient in the radial velocity across the full sky. Such a coherent pattern can only be obtained if stars are moving in the same (retrograde) direction on elongated orbits. The proper motions, depicted by the arrows, reveal a rather complex velocity field. This is expected, given the large mass of the progenitor object and the short mixing timescales in the inner Galaxy. 2 Nonetheless, in this complexity we see streams: close stars often move in the same direction. This is a very significant effect as established  \nby comparing to mock sets constructed assuming a multivariate Gaussian for the velocities (see Methods for details).  \nWe conclude that the halo near the Sun is strongly dominated by a single structure of accreted origin, as hinted also by other work, 13,14 and leaving little room for an in-situ contribution. 15 It is however, not necessarily representative of the whole stellar halo, as debris from other accreted large objects (with e.g. different chemical abundance patterns) might dominate elsewhere in the Galaxy. We also conclude that the Milky Way disk experienced a significant merger in its history. We estimate the mass-ratio of this merger at the time it took place as M GE vir M MW vir = f MW f GE × M GE ∗ M MW ∗ , where f is the ratio of the luminous-to-halo mass of the object. At the present-time, f MW, 0 ∼ 0 . 04 for the Milky Way, 29 and if we assume that Gaia-Enceladus would be similar to the Large Magellanic Cloud had it evolved in isolation, then f GE, 0 ∼ 0 . 01. 22 It has been shown 30 that the redshift evolution of f between z = 2 and z = 0 for objects of the Magellanic Cloud scale and the Milky Way is similar, implying that f MW /f GE = f MW, 0 /f GE, 0 ∼ 4. Therefore, taking M MW ∗ at the time of the merger to be the mass of the thick disk, 29 i.e. ∼ 10 10 M glyph[circledot] , we obtain a mass-ratio for the merger of ∼ 0 . 24. This implies that the merging of Gaia-Enceladus must have led to significant heating and to the formation of a thick(er) disk.", 'REFERENCES': '- 31 Schonrich, R., Binney, J., & Dehnen, W. Local kinematics and the local standard of rest. Mon. Not. R. Astron. Soc., 403 , 1829-1833 (2010)\n- 32 Gaia Collaboration, Katz, D., Antoja, T., et al. Gaia Data Release 2: Mapping the Milky Way disc kinematics. Astron. Astrophys., 616 , A11 (2018)\n- 33 Helmi, A., Veljanoski, J., Breddels, M. A., Tian, H., & Sales, L. V. A box full of chocolates: The rich structure of the nearby stellar halo revealed by Gaia and RAVE. Astron. Astrophys., 598 , A58 (2017)\n- 34 Arenou, F., Luri, X., Babusiaux, C., et al. Gaia Data Release 2: Catalogue validation. arXiv:1804.09375 (2018)\n- 35 Jean-Baptiste, I., Di Matteo, P., Haywood, M., et al. On the kinematic detection of accreted streams in the Gaia era: a cautionary tale. Astron. Astrophys., 604 , A106 (2017)  \n- 49 Deason, A. J., Belokurov, V., Koposov, S. E., & Lancaster, L. Astrophys. J. Lett., 862 , L1 (2018)\n- 50 Conroy, C., Bonaca, A., Naidu, R. P., et al. They Might Be Giants: An Efficient ColorBased Selection of Red Giant Stars. Astrophys. J. Lett., 861 , L16 (2018)\n- 51 Belokurov, V., Evans, N. W., Bell, E. F., et al. The Hercules-Aquila Cloud. Astrophys. J. Lett., 657 , L89-L92 (2007)\n- 52 Larsen, J. A., Cabanela, J. E., & Humphreys, R. M. Mapping the Asymmetric Thick Disk. II. Distance, Size, and Mass of the Hercules Thick Disk Cloud. Astron. J., 141 , 130 (2011)  \nData Availability Statement All data generated or analysed during this study are included in this published article as Source Data files.', 'Methods': "We describe here the motivation behind the name Gaia-Enceladus. In Greek mythology Enceladus was one of the Giants (Titans), and the offspring of Gaia (which represents the Earth), and Uranus (representing the Sky). Enceladus was said to be buried under Mount Etna and responsible for earthquakes in the region. The analogies to the accreted galaxy reported and characterized in this paper are many, and they include: i) being offspring of Gaia and the sky, ii) having been a 'giant' compared to other past and present satellite galaxies of the Milky Way, iii) being buried (in reality first disrupted by the Milky Way and then buried, also in the Gaia data as it were), and iv) being responsible for seismic activity  \n(i.e. shaking the Milky Way and thereby leading to the formation of its thick disk). We refer to the accreted galaxy as Gaia-Enceladus to avoid confusion with one of Saturn's moons, also named Enceladus.", 'A. Dataset, selection criteria and the effect of systematics': "For the work presented in the main section of the paper, we selected stars from the Gaia 6D-dataset 17 with small relative parallax error glyph[pi1] /σ glyph[pi1] > 5, which allows us to compute their distance as d = 1 /glyph[pi1] . For Figure 1, we consider only stars with glyph[pi1] > 0 . 4 mas (i.e. within 2.5 kpc from the Sun) to limit the impact of velocity gradients. The velocities were obtained using the appropriate matrix transformations form the observables α, δ , µ α ∗ , µ δ , v los and distances d . These velocities have then been corrected for the peculiar motion of the Sun 31 and the Local Standard of Rest velocity, assuming a value 29 of V LSR = 232 km/s.  \nWe select halo stars (such as the black points in Figure 1a) as those that satisfy | v -v LSR | > 210 km/s. 14 This condition is an attempt to remove the contribution of the disk(s), although towards the inner Galaxy, this is less effective because of the increasing velocity dispersion of disk stars. 32 To select members of the retrograde structure (such as the blue points in Figure 1a), we inspect the energy vs L z distribution of the stars in our dataset. The energy is computed assuming a Galactic potential including a thin disk, bulge and halo components. 33 For example, the left panel of Extended Data Fig. 1 shows the energy vs L z distribution for all halo stars within 5 kpc from the Sun ( glyph[pi1] > 0 . 2 mas). We have here removed stars with phot-bp-rp-excess-factor > 1 . 27 (this is enough to remove some not so well-behaved globular cluster stars so we do not apply a colour-dependent correction 34 ). This figure shows that the regions occupied by the retrograde structure and by the disk are relatively well-separated. There is however some amount of overlap, particularly for higher binding energies and lower angular momenta. Therefore even the selection criteria of L z < 150 kpc km/s and E > -1 . 8 × 10 5 km 2 /s 2 , indicated by the straightlines, will not yield a pure (thick disk free) sample of stars in the structure. This figure reveals also the large range of energies in the structure, indicating that member stars are expected over a large range of distances.  \nBecause the energies of stars depend on the gravitational potential of the Galaxy, and its form and amplitude are not so well-constrained beyond the Solar neighbourhood, we use a criterion based only on L z to find additional members of the structure/Gaia-Enceladus beyond the immediate vicinity of the Sun (as for Fig. 3 of the main section). The central panel of Extended Data Fig. 1 shows the L z vs Galactocentric distance in the disk plane R , for all stars in the Gaia 6D-dataset with glyph[pi1] /σ glyph[pi1] > 5, and including stars with parallaxes  \nglyph[pi1] > 0 . 2 mas. This plot shows that a selection based only on L z works relatively well to isolate Gaia-Enceladus stars near the Sun and also farther out in the Galaxy. However for the inner regions there is much more overlap and hence the distinction between the thick disk and Gaia-Enceladus is less straightforward, and the amount of contamination by thick disk stars is likely to be much higher. Furthermore, we expect the orbits of some stars in the progenitor of the thick disk to have been perturbed so significantly during the merger 35 that they will 'mingle' with those from Gaia-Enceladus.  \nThe rightmost panel of Extended Data Fig. 1 shows the z -angular momentum as function of cylindrical radius of stellar particles in a simulation of the merger of a pre-existing disk and a massive satellite 18,36 (the same of Fig. 1b). The example here corresponds to the redshift z = 1 simulation of a disk with M ∗ = 1 . 2 × 10 10 M glyph[circledot] and a satellite with M ∗ ,sat = 2 . 4 × 10 9 M glyph[circledot] . Because of the lower host mass used in this simulation (compared to the present-day mass of the Milky Way), the spatial scales and velocities are typically smaller compared to the data. Therefore in the simulations, we consider as solar vicinity a volume centered at R sim sun = 2 . 4 × R final thick , where R final thick = 2 . 26 kpc. 18 We also scale the positions by R sun /R sim sun = 1 . 5 and the velocities by v thick , sun /v final,sim thick , where v thick , sun = 173 km/s is the rotational velocity of the thick disk near the Sun 37 and v final,sim thick is that of the thick disk in the simulation at R sim sun . Extended Data Fig. 1c shows that like for the data, the separation between accreted and host disk stars is less effective for small radii.  \nFor Fig. 2a, we have cross-matched the catalogues from Gaia DR2 and APOGEE 9, 38 DR14 and retained only stars with estimated distances from both these catalogues (i.e. spectrophotometric and trigonometric parallaxes) consistent with each other at the 2 σ level. We also impose a relative parallax error of 20%. More than 100,000 stars within 5 kpc from the Sun satisfy these conditions. The abundances shown in Fig. 2 stem from the ASCAP pipeline. 39  \nThe presence of a parallax zero-point offset in the Gaia data 40 has been established thoroughly, and is partly (if not only) due to a degeneracy between the parallax and the basic-angle variation of the Gaia satellite. 41 Its amplitude varies with location on the sky, 34,40 and is on average -0 . 029 mas and has an RMS of ∼ 0 . 03 mas. 27 Such variations make it very difficult to perform a correction a posteriori for the full Gaia DR2 dataset (although the expectation is that its effect will be smaller for Gaia DR3). The discovery and characterization of Gaia-Enceladus was done using stars with parallaxes glyph[pi1] > 0 . 4 mas for Fig. 1 of the main section, and in Fig. 2 for stars with glyph[pi1] > 0 . 2 mas from the cross-match of Gaia and APOGEE. We therefore expect the derived kinematic and dynamical quantities for these subsets to be largely unaffected by the systematic parallax error. However, for Figs. 3 and 4 of the main section of the paper, we selected stars on the basis of their L z although we  \nfocused on properties which are independent of the parallax, such as position on the sky and proper motions. Nonetheless, to establish how important the parallax zero-point offset is on the selection via L z we perform the following test.  \nWe use the Gaia Universe Model Snapshot GUMS v18.0.0, 42 and select stars according to the following criteria: 6 ≤ G ≤ 13 . 0, 0 . 2 ≤ log g ≤ 5 and 3000 ≤ T eff ≤ 9000 K. This selection leads to a total of 7403454 stars distributed across all Galactic components. For these stars we compute error-free velocities and L z . We convolve their true parallax with a Gaussian with a dispersion of depending on the magnitude of the star. 43 The parallax is reconvolved with a Gaussian with a mean of -0 . 029 mas and a dispersion of 0.030 mas. Using these observed parallaxes, we compute 'observed' velocities and L z .  \nWe find that for measured distances smaller than 5 kpc, there is no shift in the derived L z , while for a shell between 5 and 7.5 kpc the median amplitude of the shift is ∼ -50 kpc km/s, making the observed L z more retrograde. For a shell between 9 and 10 kpc, the median shift is small and has an amplitude of 20 kpc km/s, presumably reflecting that at such large distances, the random errors on the individual stars' measurements dominate. The results are shown in the left panel of Extended Data Fig. 2 where we plot the difference between the true (initial) and 'measured' distributions of L z for stars 'observed' to be located at distances between 5 and 10 kpc, for l = ( -60 o , -20 o ). The panel on the right shows the distribution of the mean value of the difference over the whole sky, and although it reveals certain patterns, these are different from those seen in Fig. 3. As stated in the main section of the paper, the lack of distant stars in the regions outside of the contours plotted in Fig. 4, is the result of a quality cut in the relative parallax error of 20%. This selection criterion allows for parallax errors in the range 0.02 to 0.04 mas for the most distant stars (with 0 . 1 < glyph[pi1] < 0 . 2 mas), and these are only reached in those regions of the sky that have been surveyed more frequently by Gaia , such as around the ecliptic poles. The Gaia RR Lyrae stars associated to Gaia-Enceladus suffer of course also from this effect, as a lower number of visits leads to more difficult identification and hence to lower levels of completeness. 26", 'B. Random sets and significance of features': "To understand how different the dynamical properties of the Gaia 6D dataset are in comparison to a smooth distribution, we plot the distribution of velocities in Extended Data Fig. 3a and of E vs L z in Extended Data Fig. 3b for randomized datasets. These smooth datasets have been obtained from the Gaia data shown in Fig. 1a and in Extended Data Fig. 1a, respectively, by re-shuffling the velocities. That is, for each star, we assign randomly  \na v y and v z velocity from two other stars in the sample. This results in distributions with the same 1D velocity distributions as the original data, but without any correlations or lumpiness. The comparison of Fig. 1a to Extended Data Fig. 3a shows that the distribution in the random dataset is indeed much smoother than the data, and that the overall velocity dispersion in the y -direction has increased because there no longer is a clear separation between the region occupied by Gaia-Enceladus and by the thick disk. The comparison of Figs. 1b and Extended Data Fig. 3b is even more revealing and clearly shows that the structure defined in E vs L z by Gaia-Enceladus stars has effectively disappeared in the randomized dataset. Similar conclusions are reached when, instead of using a reshuffled dataset, we compare the distributions to those in the GUMS model.  \nFig. 4 of the main section of the paper shows the radial velocities and proper motions (corrected for the Solar and for the Local Standard of Rest motions) for stars with 0 . 1 < glyph[pi1] < 0 . 2 mas and -1500 < L z < 150 kpc km/s. These stars are tentative members of Gaia-Enceladus, although as discussed earlier towards the inner Galaxy contamination by thick disk stars becomes more important for large distances. The arrows depicting the proper motions suggest that stars that are closeby on the sky move in similar directions. We establish here whether this is significant by comparing to a mock dataset.  \nThe mock dataset uses the measured positions of the stars that are plotted in Fig. 4, but their velocities are generated randomly according to a multivariate Gaussian distribution with dispersions in v R , v φ and v z of 141, 78 and 94 km/s respectively. 44 During the process of generation, we only keep stars' velocities that satisfy -1500 < L z < 150 kpc km/s, as in the real data. To quantify the degree of coherence in the proper motions of neighbouring stars on the sky, we perform the following test. For each star, we find its nearest neighbour on the sky, and then determine the angle between their proper motion vectors. We then count the number of such pairs having a given angle. Extended Data Fig. 4 shows the distribution of these pairs for the Gaia subsample (in blue) and for the mock (in red). There is a clear excess of pairs of stars with similar directions of motion in the data in comparison to the mock.", 'C. Context and link to other substructures': "Hints of the presence of a population like Gaia-Enceladus have been reported in the literature in the last two decades, and were typically based on small samples of stars. These hints were of chemo-dynamical nature 6, 8, 45 and sometimes attributed to accretion, 4, 46 but also based purely on chemical signatures, such as the α -poor sequence. 7, 19 More recently, cross-matches to the first data release of the Gaia mission 47 also revealed the contrast between  \nthe metal-rich population supported by prograde rotation and associated to the tail of the thick disk, 48 and the metal-poor halo, i.e. what we have just identified as Gaia-Enceladus. Furthermore, in one study 13 the difference in the kinematics of these two populations, and the measurement of a very radially biased velocity ellipsoid for halo stars with [Fe/H] > -1 . 7, led to the proposal that this population (which was termed 'Gaia sausage') could be the result of a significant merger. Although this could be also attributed to an in-situ formation via a radial collapse, this hypothesis gained further supported by their orbits leading to the break in the halo density profile at ∼ 20 kpc. 49 All of these pieces together outline the case for the discovery and detailed characterization of Gaia-Enceladus reported here.  \nThe more distant Gaia-Enceladus debris occupies large portions of the sky not extensively covered by other existing surveys. There is however, a recent detection of an overdensity identified in PanSTARRS and WISE with the help of Gaia proper motions, 50 which overlaps with the northern part of the more distant Gaia-Enceladus stars for -2 < µ α < -1 mas/yr and -1 < µ δ < 0 mas/yr, and partly (but not fully because of the PanSTARRS footprint) with the southern part, for 0 < µ α < 1 mas/yr and -3 < µ δ < -1 mas/yr. There could potentially be also a relation to the Hercules Aquila Cloud 51 identified in SDSS, although this appears to be offset both in the northern and southern hemispheres and located at a larger distance. The location on the sky of intermediate distance Gaia-Enceladus stars would seem to overlap with the Hercules thick disk cloud, 52 especially in the fourth Galactic quadrant below the Galactic plane.", 'Extended Data Figure Captions': "Extended Data Figure 1. Slices of phase-space used to isolate Gaia-Enceladus stars . Panel a): Energy E vs L z for stars in the 6D Gaia dataset, satisfying the quality criteria described in the text, with glyph[pi1] > 0 . 2 mas (5 kpc from the Sun) and with | v -v LSR | > 210 km/s. The straightlines indicate the criteria used to select Gaia-Enceladus stars, namely -1500 < L z < 150 kpc km/s and E > -1 . 8 × 10 5 km 2 /s 2 . These criteria follow roughly the structure's shape (see for comparison Extended Data Fig. 3b), but are slightly conservative for the upper limit of L z to prevent too much contamination by the thick disk. However, small shifts such as considering an upper limit of 250 kpc km/s or a lower limit of -750 kpc km/s for L z , or E > -2 × 10 5 km 2 /s 2 do not result in drastic changes to the results presented in the paper. Panel b): L z vs Galactocentric distance R for all stars in the 6D Gaia with glyph[pi1] > 0 . 2 mas. The black points are the halo sample shown in panel a). Panel c): same as panel b) for star particles in the merger simulation 18 shown in Fig. 1b, where blue correspond to the stars from the satellite, and grey to the host disk, and the positions and velocities have been scaled as described in the text. In this figure, the energy has been scaled by E sun (which is -1 . 63 × 10 5 km 2 /s 2 in the Galactic potential used), L z by L z, sun = 1902 . 4 kpc km/s, and R by the solar distance R sun = 8 . 2 kpc.  \nExtended Data Figure 2. Effect of a zero-point offset in the parallax on L z . Panel a) shows the distribution of the difference between the initial and 'measured' (after error convolution) L z for GUMS stars with 'measured' distances between 5 and 10 kpc and with l = ( -60 o , -20 o ). Panel b) shows the mean value of the difference over the full sky.  \nExtended Data Figure 3. Distribution of stars' dynamical properties for a smooth dataset . Panel a) shows the velocity distribution and panel b) the E vs L z distribution for a dataset obtained by reshuffling the velocities of the stars plotted in Fig. 1a and in Extended Data Fig. 1a, respectively. The visual comparison to those figures shows that these random sets are less clumped than the observed distributions of the Gaia halo stars.  \nExtended Data Figure 4. Distribution of angles between the proper motion vectors for neighbouring stars on the sky. The blue and red histograms correspond respectively, to Gaia-Enceladus and to a mock dataset. This mock dataset uses the positions of the stars in Gaia-Enceladus, but velocities generated randomly according to a multivariate Gaussian distribution, 44 after which only stars' velocities that satisfy -1500 < L z < 150 kpc km/s are kept, as in the real data. For each star, we find its nearest neighbour on the sky, and then determine the angle ∆ θ between their proper motion vectors for the data and for the mock. We then count the number of such pairs having a given angle ∆ θ .", 'Extended Data Figures': 'Extended Data Fig. 1  \n<!-- image -->  \nExtended Data Fig. 2  \n<!-- image -->  \n<!-- image -->  \nExtended Data Fig. 3  \n<!-- image -->  \nExtended Data Fig. 4  \n<!-- image -->'}