Streams, substructures and the early history of the Milky Way

The advent of Gaia's 2nd data release in combination with large spectroscopic surveys are revolutionizing our understanding of the Galaxy. Thanks to these and the knowledge accumulated thus far, a more mature picture of the evolution of the early Milky Way is emerging: * Two of the traditional Galactic components, i.e. the stellar halo and the thick disk, appear to be intimately linked: stars with halo-like kinematics originate in similar proportions, from a"heated"(thick) disk and from debris from a system named Gaia-Enceladus. Gaia-Enceladus was the last big merger event experienced by the Milky Way and probably completed around 10 Gyr ago. The puffed-up stars now present in the halo as a consequence of the merger have thus exposed the existence of a disk component at z ~ 1.8. * The Helmi streams, Sequoia, and Thamnos are amongst the newly uncovered or better characterized merger events. Knowledge of their progenitor's properties, star formation and chemical histories is still incomplete. * Debris' from different objects often overlap in phase-space. A task for the next years will be to use spectroscopic surveys for chemical labelling and to disentangle events from one another using dimensions other than only phase-space, metallicity or [alpha/Fe]. * These surveys will also provide line-of-sight velocities missing for faint stars and more accurate distance determinations for distant objects. The resulting samples of stars will cover a much wider volume of the Galaxy allowing, for example, linking kinematic substructures in the inner halo to spatial overdensities in the outer halo. * All the results obtained so far are in-line with expectations of current cosmological models. Yet, tailored hydrodynamical simulations as well as"constrained"cosmological simulations are needed to push our knowledge of the assembly of the Milky Way back to the earliest times. [abridged]


INTRODUCTION
The current are very exciting times for research on streams and substructures, and their use to shed light onto the early history of our own galaxy, the Milky Way. Although the field now known as Galactic Archaeology has a long history, it is hard to overstate the impact of the second data release from the Gaia Mission (DR2, Gaia Collaboration et al. 2018b), which took place on April 25th 2018. The combination with data already available from many large spectroscopic surveys such as APOGEE 1 (Majewski et al. 2017), GALAH 2 (De Silva et al. 2015), RAVE 3 (Kunder et al. 2017), and LAMOST 4 (Deng et al. 2012), has helped to obtain a much clearer picture of how the Milky Way, and in particular its older components, have evolved since z ∼ 2, or equivalently 10 Gyr ago.
These new datasets are allowing putting together and in a broader context, the many pieces of the puzzle previously reported in the literature to give a much more complete view of the Galaxy's past. The current generation is quite fortunate to be part of this chapter in the history of Galactic astronomy. It is very exciting that we might actually know how and when the Milky Way experienced its last big merger and that it seems likely that this event gave rise to most of the halo near the Sun, which would be predominantly be composed of debris from a single object that was accreted about 10 Gyr ago, and of heated disk present at the time. This is what Gaia has unraveled in conjunction with high resolution spectroscopic surveys, particularly with APOGEE. The rapid progress made in the field since Gaia DR2 has been possible thanks to the work of many scientists also before DR2, as their work allowed deriving relatively quickly a rather clear, although not yet fully settled, picture of the sequence of events. This is in fact an example of one of the pillars of the scientific enterprise: that we build on previous knowledge. It would have taken much longer to pin down Galactic history to the extent reached thus far had these earlier works not been carried out. The 2nd data release of the Gaia mission, even if only based on data taken during less than half of the mission's nominal lifetime (22 months out of 60), has really helped us to move from a fragmented view to seeing Galactic history in full glory.
Many excellent reviews have been written over the past 20 years on Galactic archaeology and near-field cosmology starting with Freeman & Bland-Hawthorn (2002); and include also Frebel & Norris (2015) on first stars and their use for (near-field) cosmology; on the structure and dynamics of the Galaxy by Bland-Hawthorn & Gerhard (2016); on substructure and tidal debris by Belokurov (2013) and Johnston (2016), as well as the introduction to the Galactic halo by Helmi (2008). An interesting exercise is to read the reviews using the information that we have recently acquired about our Galaxy. The reader is encouraged to put on the new Gaia glasses when going through the findings reported in those studies. Hopefully s/he will note that there is much consistency in the results obtained so far, and hopefully also these reviews will aid the readers in constructing their own narrative on the basis of the information and hints that we had but we did not fully understand at the time.
The first objective of this review is thus to present the state-of-the-art in the context of what was previously known about our Galaxy. It should be noted that because we are still in the process of digesting the most recent results from the many ongoing surveys focused on the Milky Way, and because there is much more data to come in the next 5 to 10 years, it is particularly challenging to give an overview that is complete and that will stand the test of time. The emphasis and sometimes the interpretation of the recent discoveries reflects this author's own perspective and understanding, while still aiming for an objective and solid account of the facts.
Another objective of this review is also to bring out new venues for research now that we have a much better, albeit also sketchy and as just acknowledged, still in a state of flux understanding of the assembly of the Milky Way. As described especially in the second half of this review, there are still many small and not so small details missing. Solving these will require substantial effort. We will need more detailed modeling and better hydrodynamical and cosmological simulations. We will have to assemble large high resolution spectroscopic datasets with the chemical abundances of millions of stars to be able to pin down their site of formation, to be able to label as it were, the stars' origin. It should be possible to go back in time even further than 10 Gyr ago, perhaps out to redshift 6 -10 by studying stars in the different structures of the Milky Way.
This review starts in Sec. 2 with a brief description of the different Galactic components following a traditional approach. In Sec. 3 we move on to Galactic archaeology, and discuss the fossils and tools that are available to do this type of work. Then we dive in Sec. 4 into one of the components that holds clues to the evolution of the Galaxy at early times, namely the Galactic stellar halo. We describe the most recent discoveries and how they link to the formation of another ancient component, the thick disk. We focus on this latter component in Sec. 5. In this journey, we describe not just the data but we also discuss predictions from simulations and models. In Sec. 6 we describe the next steps, those that would seem to be necessary to really fully unravel how the Milky Way was put together.
These as well as the most important conclusions are summarized in Sec. 7.

Brief description
The Milky Way is, in general terms, a fairly typical disk galaxy (Bland-Hawthorn & Gerhard 2016). Its estimated stellar mass is ∼ 5 × 10 10 M , which implies a luminosity close to the characteristic value L * of the galaxy luminosity function. Given its circular velocity of Vmax ∼ 240 km/s (see e.g. Gravity Collaboration et al. 2019), it may be slightly subluminous as it lies a bit below -but within 1σ, the Tully-Fisher relation.
The Milky Way has several visible components: a thin disk, thick disk, bulge/bar and a stellar halo as shown in Figure 1. Each of these components has individual characteristics. Not only do their stars differ in their spatial distribution but of course, also kinematically as shown in Figure 2. Furthermore their age and chemical distributions are also different. This implies that the components are truly physically distinct. Their constituent stars inform us about the various processes that are important in the build up of a galaxy throughout its life. The Milky Way and its various components. This image was obtained using data from the 2nd data release of the Gaia mission (Gaia Collaboration et al. 2018b). Credits: ESA/Gaia/DPAC, CC BY-SA 3.0 IGO.
2.1.1. Short summary of the main characteristics of the Galactic components.
• The thin disk is the site of ongoing star formation and it is the most characteristic component of the Galaxy, giving it its name. Its current star formation rate is es- Velocity distribution of stars in the solar neighborhood as determined by Gaia. In this figure, all stars from Gaia DR2 with full phase-space information, located within 1 kpc from the Sun, and with relatively accurate parallaxes, i.e. with /σ ≥ 5 have been considered. The nearby halo stars are plotted with black dots and defined as those that satisfy |V − V LSR | > 210 km/s, for V LSR = 232 km/s. The blue density maps reveal the contribution of the thin and thick disks. The "banana"-shaped structure seen in the left panel reveals an important contribution of "hot" thick disk-like stars to the halo. Credits: H.H. Koppelman (see also Fig. 2 in Koppelman et al. 2018).
timated to be ∼ 1.6 M /yr (Licquia & Newman 2015), and it seems to have been forming stars at least for 8 or 9 Gyr (Tononi et al. 2019). It is rotationally supported, and most stars move on fairly circular orbits. • The thick disk is a thicker, more diffuse and hotter component than the thin disk.
Its stars are older than the oldest stars in the thin disk, with estimates using white dwarfs in the Solar vicinity suggesting by at least ∼ 1.6 Gyr (Kilic et al. 2017). Its metallicity distribution function peaks at a lower metallicity value of [Fe/H] ∼ −0.5, and its stars define a separate chemical sequence in e.g. [α/Fe] vs [Fe/H] space, from that defined by the thin disk (Bensby et al. 2003, Fuhrmann 2011, which can be attributed to a different (shorter and more intense) SFH (see e.g. Chiappini et al. 1997. We discuss it in more detail in Sec. 5. • The bar/bulge is the most centrally concentrated component, and because it is heavily obscured, our current understanding is somewhat limited, although significant progress has recently been made thanks to new surveys, as described in e.g. the reviews by Barbuy et al. (2018) and Zoccali (2019). The presence of a "classical" bulge (i.e. of spherical shape, formed quickly, dispersion supported) is still debated but its contribution has been constrained by the observed kinematics to be small (< 8% of the mass of the disk, Shen et al. 2010). Most of the bulge is in a rotating triaxial structure, the Galactic bar, whose orientation, pattern speed and in particular, exact extent have undergone revision lately, where recent work suggests a rather long bar . Spectroscopic studies show a mix of populations present in the central regions (Ness et al. 2013), some of which are very old (more than 13 Gyr), and metal-rich with [Fe/H] values up to +0.5 dex, and some that resemble other Galactic components, such as the thick disk and stellar halo, all of which of course, peak in terms of their spatial density in the inner Galaxy. • The stellar halo is the most extended component but at the same time it is rather centrally concentrated: the half-light radius traced by the metal-poor globular clusters is ∼ 0.5 kpc (Bica et al. 2006). It is oblate in the inner regions with q ∼ 0.6, and its density is well-modelled by a broken power-law (Deason et al. 2011, Xue et al. 2015.
The most recent estimates of its total mass yield ∼ 1.3 × 10 9 M , Mackereth & Bovy 2020. The stellar halo contains very metal-poor and old stars. It will be discussed in detail in Sec. 4 of this review. • The above items refer to the stellar components of the Galaxy, but there is also warmionized gas (in a halo or circum-galactic medium, Richter 2017, Zheng et al. 2019, and cold gas mostly in the disk.
If our understanding of Gravity is correct, the Galaxy is embedded in a dark matter halo, where most of the mass of the system is located. The characteristics of this halo are not very well constrained. Current estimates of its mass based on Gaia DR2 by Posti & Helmi (2019), Watkins et al. (2019) give ∼ 1.3 × 10 12 M (consistent with the range of values quoted in Bland-Hawthorn & Gerhard 2016). Its shape is uncertain and has been the subject of significant debate (Ibata et al. 2001, Helmi 2004, Johnston et al. 2005, Ibata et al. 2013. It is likely slightly oblate in the central regions (Koposov et al. 2010, although Wegg et al. (2019 argues for spherical) and changes to a triaxial shape at large distances (Law & Majewski 2010, Vera-Ciro & Helmi 2013, with the longest axis in the direction perpendicular to the disk (Banerjee & Jog 2011, Vera-Ciro & Helmi 2013, Bowden et al. 2016, Posti & Helmi 2019. The density profile of the dark halo has received less attention thus far (but see Taylor et al. 2016, Fardal et al. 2019, Eadie & Jurić 2019, Yang et al. 2019). An interesting question is the degree of lumpiness of the mass distribution and whether it is consistent with expectations from cold dark matter simulations, which predict myriads of (dark) satellites (Moore et al. 1999, Klypin et al. 1999. Recent work on streams is beginning to reveal a complexity that may require the consideration of perturbations by e.g.

Link between the components and physical processes in galaxy evolution
The differing characteristics of the various Galactic components suggest each had its own formation path. Nonetheless it is likely that these paths were interlinked. It should largely be possible to unravel these formation paths using stars, since these retain memory of their origin. This idea constitutes the pillar of Galactic archaeology, as we discuss in greater detail in Sec. 3.
The Λ-cold dark matter (ΛCDM) model provides a framework to understand how galaxies form and evolve from first principles (see e.g. the review by Frenk & White 2012). In this model, galaxies form inside dark matter halos (White & Rees 1978). Most of the dark halos properties (such as their mass function, abundance number, et cetera) depend on characteristics of the cosmological model, including for example the power spectrum of density fluctuations, the type of dark matter, and the values of the cosmological parameters (as discussed extensively in the book by Mo et al. 2010). Because in the concordance model there is ∼ 6× more mass in dark matter than in baryons (this is supported by measurements of e.g. the fluctuations in the cosmic microwave background, Planck Collaboration et al. 2016), many of the properties of galaxies such as how they cluster or their dynamics, are largely dictated by their dark halos. A direct example of this is the process of halo collapse and formation, during which dark halos attract baryonic material from which the (visible components of) galaxies can form. For the baryons in a gaseous configuration to be able to cool and form stars, several conditions need to satisfied (dictated by e.g. cooling and heating processes, and dynamical timescales, see Mo et al. 2010). If these conditions are satisfied, the gas will cool and collapse to the center of their halos while conserving some amount of angular momentum. This results in a gaseous disk that is rotationally supported (Mo et al. 1998), with some amount of random motion depending on the state of the gas (Bournaud et al. 2009). Note that particularly in the early universe, gas can also be directly accreted as a cold flow and feed the forming galaxy (Dekel et al. 2009). In the cold gas disk, stars will start to form. In fact, most star formation in the Universe takes place quiescently and is not associated to large starbursts (see Brinchmann et al. 2004, Elbaz et al. 2011. In the ΛCDM model structure formation proceeds hierarchically, via mergers. At early times mergers were more frequent because of the higher density of the Universe. This means that galaxies were more prone to merge with other galaxies, and hence their disks were more vulnerable. Depending on the mass ratio, such an event could lead to the formation of a bulge (Barnes 1992), or merely to the thickening of the disk (Quinn et al. 1993), and possibly also to the formation of a halo of stars from the original disk and from the destroyed satellite (Zolotov et al. 2009, Purcell et al. 2010, as seems to have happened for the Milky Way (see Sec. 4.2 for details). Depending on the characteristics of the merger, such an event could have also triggered the formation of a bar (Gerin et al. 1990). It is in fact likely that the Galactic bar has originated from a disk instability. However it is not clear whether the bar had its origin in the thin disk (Martinez-Valpuesta & Gerhard 2013), or whether the metallicity gradient seen in the bar implies that some of the stars have their origin in the thick disk (see Di , Fragkoudi et al. 2018, and references therein), as suggested also by their similar chemical abundance patterns (Alves-Brito et al. 2010).
These examples show that there may be strong links between different components of the Milky Way, and that some could even have gotten their current configuration during or triggered by the same event. Furthermore these components may share a fraction of their stellar populations, such as the bar with the (primordial) thick disk, or the halo with the primordial thick disk. On the other hand galaxies at earlier times had higher gas fractions, which could also imply that mergers may have indirectly led to the formation of a significant stellar population in a Galactic component via the triggering of a starburst, as perhaps was the case for the thick disk (see Gallart et al. 2019, and Sec. 5.2 for more details).
These considerations highlight why we should probably not think of our Galaxy in terms of separate and independent components that have no connection to each other. Rather we should aim to establish if and how they may be related given our ultimate goal to unravel the sequence of events that took place in the history of the Milky Way.

Introduction
Today's commonly used phrase Galactic archaeology is often applied to describe research on the formation and history of the Milky Way and its stellar populations. The work of Roman (1950, and several of her subsequent papers) showing that stars with different chemistry also have different kinematics, has been recognized to have been very influential 5 . The papers by Eggen et al. (1962) and Searle & Zinn (1978), as well as Tinsley (1980) more generally for galaxies, can arguably be considered as pioneering the field.
In its modern form the idea behind "Galactic archaeology" is to use the properties of long-lived stars to reconstruct the history much in the same way archaeologists use artifacts or "rubble", to learn about the past. Possibly one of the first printed records of the use of the word "archaeology" in an astronomical context is an article in The ESO Messenger by Spite & Spite (1979), where there is a reference to "astro-archaeology". In this paper the authors aim to use old stars to understand the build-up of metals in the universe. The term "Galactic archaeology" in a more dynamical context is used in the Report of IAU Commission 33: Structure and dynamics of the Galactic system (Burton 1988), in Section 13 (in charge of J. Binney): "Perhaps it is not too fanciful to imagine a field of galactic archaeology opening up, in which painstaking sifting of the contents of each element of phase-space will enable us to piece together a fairly complete picture of how our Galaxy grew to its present grandeur and prosperity".
The turn of the century is approximately the time that the phrase "Galactic archaeology" was adopted widely by the community, as it begins to appear more frequently both in talks as in the printed literature, also in part because of the very influential reviews by Bland-Hawthorn & Freeman (2000), Freeman & Bland-Hawthorn (2002) (see also Bland-Hawthorn 1999, who introduced the term "near-field cosmology"). Impetus to the field was undoubtedly given by the discovery by Ibata et al. (1994) of the Sagittarius dwarf as direct evidence of an ongoing merger, and to some extent subsequently by discovery of debris streams near the Sun from a past merger in the Hipparcos 6 data (Perryman et al. 1997) by .
This time also coincides with the maturing of galaxy formation models (Kauffmann et al. 1993, Baugh et al. 1998, Somerville & Primack 1999, and the establishment of the ΛCDM model as the concordance cosmological model. This allowed significant progress in the theoretical predictions concerning what a galaxy like the Milky Way should have experienced in its lifetime. Thus "Galactic archaeology" could also be guided by theory and some aspects of the cosmological models could now tested directly from the perspective of the Milky Way. This spirit is particularly evident in the 3rd Stromlo symposium on the Galactic halo edited by Gibson et al. (1999) which took place in Canberra in 1998. For example, the article describing the conference highlights (de Zeeuw & Norris 1999), as well as a quick inspection of the index of the proceedings will reveal that the theme of accretion and mergers and the use of the fossil record to reconstruct Galactic history was present in many of the participants' contributions to the meeting.
What does the "Galactic archaeology" actually mean? As already mentioned the idea behind it is that stars have memory of their origin. Low-mass stars live longer than the age of the universe, and hence some will have formed at very early times and have survived until the present day. They will have retained in their atmospheres a fossil record of the environment in which they were born. This is because the chemical composition of a star's atmosphere, particularly if it has not yet evolved off the main sequence, reflects the chemical composition of the interstellar medium (the molecular cloud) in which it formed. This means access to the physical conditions present at the time of formation of the star. For very old stars, the conditions might have been very different than today (leading for example to different initial mass functions, etc), and therefore such stars provide us with a window into the early Universe (Frebel & Norris 2015). On the other hand, stars with similar chemical abundance patterns likely have a common origin. The search for common "DNA" would then allow the identification of stars with a similar history, and is known as chemical "tagging". The foundations of this approach were put forward by Freeman & Bland-Hawthorn (2002) and are briefly discussed in Sec. 3.2.
Another particularly useful way to track Galactic history is through precise measurements of stellar ages. Knowledge of the ages of stars would permit dating the sequence of events that led to the formation of the different components of the Galaxy. However obtaining precise ages for very old stars is very difficult. Even 10% errors at 10 Gyr imply going from redshift 1.8 to 2.3, and a difference of only 2 Gyr exists for a star born at redshift 2 or at redshift 6. Nonetheless the combination of ages and chemical abundances of stars is very powerful and can be used to establish a timeline (i.e. in a closed system, stars born later will be more metal-rich).
Stars also retain memory of their origin in the way they move. For example as a galaxy gets torn apart by the tidal forces of a larger system like the Milky Way, the stripped stars continue to follow similar trajectories as their progenitor system (Johnston et al. 1996, Johnston 1998). This implies that if the Milky Way halo is the result of the mergers of many different objects, their stars should define streams that crisscross the whole Galaxy ). As will become clear in Sec. 3.3, access to full phase-space information is critical to reconstruct the past history of the Galaxy using dynamics.
3.2. Astrophysical properties of stars: Chemical abundances and ages as a tool 3.2.1. Chemical abundances. The discovery that stars with different metallicity (or iron abundance) have different chemical abundance patterns was first hinted at the end of the 1960s (Conti et al. 1967), and one of the first systematic studies of metal-poor stars is the work of Sneden et al. (1979), and interpreted in the context of supernovae type I and II and Galactic nucleosynthesis models.
The reason for the variety of chemical elemental abundance patterns is that different elements are produced in different environments and on a range of timescales (McWilliam 1997). For example, α-elements such as O, Mg, Si, Ca, S, and Ti are released in large amounts during the explosion of a massive star as a Supernova, an event that occurs only a few million years after the star's birth. On the other hand, iron-peak elements are also produced in supernovae (the so-called Type Ia) that are the result of a thermonuclear explosion of a white dwarf in a binary system, although the details of the burning and the number of white dwarfs involved or their masses are under debate (see e.g. the review Maoz et al. 2014). Because both stars in the binary are of lower mass, these SN explosions take place typically on a longer timescale than for SNII, of the order of 0.1 to a few Gyr (e.g. Matteucci & Recchi 2001). In terms of the chemical evolution of a (closed) system, we thus expect that [α/Fe] will eventually decrease as time goes by and the ISM of the system becomes polluted by SNIa. When a significant number of such explosions have occurred, the initial nearly constant [α/Fe] trend with [Fe/H] bends over and this leads to the appearance of a "knee" after which [α/Fe] can only decrease further (unless there is some fresh gas infall).
Heavier elements beyond the iron-peak are created by neutron capture processes, through the so-called slow (s) and rapid (r) processes. When the neutron flux is relatively low, i.e. the timescale between neutron captures is large compared to that of the β-decay, the s-process can occur. This can take place for example in the envelopes of Asymptotic Giant Branch (AGB) stars (Busso et al. 1999), and the contribution of low mass AGB stars (of 1 to 3 M ) appears to be particularly important in the chemical history of the Galaxy (see e.g. Bisterzo et al. 2010, and references therein, and also Battistini & Bensby 2016). For low metallicities however, stars with such masses have not had enough time to reach the AGB phase to be significant contributors of these elements (see Travaglio et al. 2004, and for a comprehensive review on s-process elements, see also Käppeler et al. 2011). A prime example of an element for which the s-process is dominant at [Fe/H] −1.5 is Ba (Arlandini et al. 1999).
The r-process, on the other hand, occurs when the neutron flux is sufficiently high to allow for rapid neutron captures. This could occur in SN II environments for example, but also in the mergers of two neutron stars and of neutron stars with black holes, or in magneto-rotational supernovae (as explored for example in the galactic simulations of Haynes & Kobayashi 2019). The recent discovery of strontium in the spectra of the kilonova following the gravitational-wave event GW170817 (Watson et al. 2019), clearly demonstrates that r-process does occur in neutron star mergers. Nonetheless, the exact sites and conditions under which the various neutron-capture elements are produced, particularly at very low metallicities, have not yet been fully settled and there may be different channels for producing them (see the excellent review of Sneden et al. 2008, and the more recent extensive review by Cowan et al. 2019). Besides strontium, a very typical r-process element is Eu, while for example Nd is produced almost equally by r-and s-processes at the solar metallicity (Arlandini et al. 1999).
The information that can be obtained from detailed chemical abundances analysis is what underpins the principle of chemical tagging, as put forward by Freeman & Bland-Hawthorn (2002). The chemical DNA of stars born in a variety of environments, will be different (De Silva et al. 2015). Although in principle each molecular cloud will have its own chemical composition, and this is likely to differ from cloud to cloud in a galaxy, in practice the differences for clouds collapsing at the present day may be small, making it very difficult to disentangle (relatively young) groups of stars of common origin on the basis of their chemistry only unless extremely accurate measurements of many different elements are available (although not impossible, see e.g. De Silva et al. 2006). It would be very interesting to associate each star in a galaxy to its parent molecular cloud because this would potentially reveal the physical processes acting on 1 -100 pc scales, i.e. the regime of the interplay between dynamics, star formation and stellar feedback. Yet this is very challenging, and thus a less demanding form of chemical tagging, known as weak tagging or chemical labelling 7 has been put forward. With "chemical labelling" we study different (larger) regions (or components) in the Galaxy to unravel for example migration mechanisms in the disk(s) (Minchev et al. 2017, Ness et al. 2019). This allows establishing whether a star now part of the thick disk actually formed in the thin disk in the inner Galaxy and migrated to the Solar vicinity (Schönrich & Binney 2009). Chemical abundances of stars in four dwarf spheroidal galaxies (in color) and in the Milky Way (in black). The left panels show the behavior with [Fe/H] of two α-elements, namely Ca and Mg, while the panels on the right correspond to the trends followed for Ba (an s-process element at [Fe/H] −1.5 or thereabouts) and Eu (r-process) with metallicity. Credits: Figure  .
In what pertains to the halo, the underlying thought behind chemical labelling is that stars born in different systems (accreted galaxies or in the proto-Milky Way) will follow their own distinct chemical sequences, because each system had its own particular star formation and chemical enrichment history. This is in fact what we see for stars associated to the different dwarf galaxies in the Local Group, as shown in Figure 3 from Tolstoy et al. (2009). Notice also how the sequences followed by the stars in the different galaxies appear to be sorted according to the mass of the system. In particular, the trend of [α/Fe] with [Fe/H] could be an interesting discriminator of stars born in accreted dwarf galaxies. Low mass galaxies that have only formed one generation of stars will likely only have high [α/Fe] at low [Fe/H], while those galaxies that have managed to sustain star formation longer, might have very low [α/Fe] even at low [Fe/H] because of inherently inefficient star formation, and hence their debris may be more easily identifiable.
Other potentially promising chemical labels for the identification of stars born in accreted dwarf galaxies appear to be r-process element abundances (see e.g. Xing et al. 2019). In the Galactic halo, there is a large scatter in [r-process/Fe] at low metallicity (as seen to some extent in the bottom right panel of Figure 3), which could be indicating a range of birth places. While most ultra-faint dwarf galaxies appear to be deficient in r-process elements, the Reticulum II galaxy, contains in proportion many r-process enhanced stars (∼ 78% compared to less than 5% in the Galactic halo, Ji et al. 2016). The way to understand this is that the events leading to the formation of r-process elements are so rare, that they have not occurred in most ultra-faint dwarfs (given their low mass). But if one event does happen, it immediately enriches the entire galaxy. Thus stars with extreme r-process abundances could have their origin in such galaxies (Brauer et al. 2019). For more massive galaxies, clustering in r-process elemental abundances might be expected (Tsujimoto et al. 2017), which combined with the behaviour of [α/Fe] or [Fe/H], could enhance their utility for chemical labelling (Skúladóttir et al. 2019).

Ages.
In comparison to chemical abundance estimation, the determination of precise ages, particularly for old stars, is much more difficult. Age determination has traditionally been done via isochrone fitting. Recently Bayesian inference tools have been employed to derive ages for large numbers of stars by using not only multi-color photometry, but also astrometric data from Gaia and chemical abundance information provided by large spectroscopic surveys (see for example Sanders & Das 2018, Queiroz et al. 2018, and also Mints & Hekker 2018. Such ages tend to be more reliable, particularly in comparison to those based only on photometry. Recently a new way of estimating ages using information about the internal structure of a star (other than the Sun) has become possible via asteroseismology. This field is growing at a fast rate and providing new insights and understanding on stellar evolution, and as a consequence also on age determination (Michel et al. 2008, Chaplin et al. 2014). Asteroseismology uses time series of photometry of outstanding quality with campaigns that may take several months or years depending on the type of star (main sequence, RGB or AGB). The photometric variations are due to internal oscillations, and their frequencies depend on the star's mass, radius and effective temperature. Because the frequencies relate in different ways to each of these parameters, the mass of a star can in principle be derived with knowledge of the basic frequencies as well as of its temperature from, for example, broad-band photometry. The knowledge of the star's mass can then be used to determine its age using stellar evolution models .
In the recent past COROT 8 (Auvergne et al. 2009) and Kepler 9 (Gilliland et al. 2010) have been providing new "golden rulers or standards" that allow for better age determination from the frequencies of oscillations of the stars. By calibrating on these, it is possible to obtain independent constraints on e.g. the gravity of a star (log g), which can then be used as a prior for the analysis of spectroscopic surveys. This then results in a larger sample of stars that have been (indirectly) calibrated, and translates into more accurate stellar parameters determinations, which in combination with isochrone fitting can then yield ages for large samples of stars (see e.g. Valentini et al. 2017). The recently launched TESS 10 (Ricker et al. 2015) and the upcoming PLATO 11 mission (Rauer et al. 2016) will monitor and characterize large samples of stars for which ages will then be readily available, plausibly much more accurately than has ever been possible until now as argued by Kollmeier et al. (2019).

Kinematical properties of stars: Dynamics as a tool
As mentioned earlier, when a galaxy is disrupted by tidal forces, its stars continue to follow closely the trajectory of the system they used to belong to. A regular orbit (a trajectory), may be characterized by the integrals of motion, such as energy E, total angular momentum (for a spherical system) or one of its components (in the case of an axisymmetric galaxy, Lz), or by the associated actions, such as JR, J φ and Jz for an axisymmetric system (Binney & Tremaine 2008). Since a small galaxy may be seen as an ensemble of stars with similar positions and velocities, this implies that also their integrals of motion (or their orbits) are similar. Hence if these are conserved through time (as expected to hold to first order for a collisionless system such as a galaxy), this implies that the tidally stripped stars will follow very similar orbits as their progenitor. This then results in the formation of a stream . A stream may thus be seen as a portion of an orbit populated by stars (to first order, see Sanders & Binney 2013, for the caveats). This explains why streams are long and narrow if originating in a small system or if formed recently (see the excellent review by Johnston 2016, for more information).
In the case of a more massive object, tides act in the same way, but the stars that are stripped at any given point in time have a larger range of values of the integrals (e.g. of energies), which results in a broader population of orbits, and hence in broader streams (and sometimes even hard to distinguish spatially). The process is not different, but the end-product has a different visual appearance and higher complexity. Furthermore, the morphological features of the debris depend also on the type of orbit of the system (Hendel & Johnston 2015, Amorisco 2017. For example, if the orbit of the progenitor is fairly radial, then shells are very pronounced. These correspond to the turning points of the orbits of the stars , Tremaine 1999. If the orbit is circular, then there are no turning points, and hence no shells. An example of the spatial evolution of debris on a somewhat radial orbit (with apocenter/pericenter ≈ 4.5) is shown in the top panels of Figure 4.
The properties of a stream depend on the extent of the parent object, the time since it formed (i.e. since a star became unbound) t and the characteristic orbital timescales, which we denote as t orb . The density of a stream at a given point in space may roughly be expressed as ρ ∝ (t orb /t) 3 × 1/(Rσ 2 ) (see , for details and the full derivation). Here R and σ are the characteristic size and velocity dispersion of the progenitor system. The dependence on time t is related to the form of the potential and the number of independent orbital frequencies (see Vogelsberger et al. 2008). Here it is assumed to be axisymmetric and the orbit to be quasi-regular (and non-resonant), hence the dependence on t −3 . The expression shows that in the first stages of the dispersal (t ∼ t orb ), the debris will have a high density and therefore remain spatially coherent, leading to easily detectable overdensities on the sky, such as those discovered in the Sloan Digital Sky Survey (SDSS 12 ) by Belokurov et al. (2006). This is typically the regime of streams orbiting in the outer halo, since there t orb is large and the tidal forces are less strong, implying also that t is small. For the inner halo, however, the orbital timescales are short, and therefore the density will decrease quickly, also for streams originating in small objects.
The behavior of stars in a stream is different if their orbits are irregular or chaotic. In that case the rate of divergence will no longer be a power-law but exponential, and phase-mixing is therefore much faster, see Price-Whelan et al. (2016). On the other hand, if the orbit is resonant stars take longer to spread out and the debris can remain spatially coherent over more extended timescales.
As time goes by, debris streams mix spatially, i.e. they become long enough that they may cross each other, and therefore one single system can be responsible for multiple streams in a given location in the Galaxy. What characterizes each of the streams is that locally the stars have very similar velocities (in this sense, the stars are truly streaming through the host galaxy). Furthermore, because of the conservation of phase-space density (or volume), as a stream becomes longer and longer, its velocity dispersion will have to decrease, i.e.

Figure 4
Comparison of various spaces commonly used to identify merger debris, namely (a) frequency space; (b) velocity space; (c) energy and Lz; (d) orbital pericenter vs apocenter; (e) actions space and (f) phase-space slice of r vs vr. The accreted satellite depicted here was evolved in a spherical Plummer potential (of mass 10 12 M and b ∼ 22 kpc) for 10 Gyr. It was non-self-gravitating, spherical and represented with a 6D-Gaussian with σx ∼ 1 kpc and σv = 22 km/s. These characteristics make it comparable to the dwarf elliptical galaxy NGC185 and whose luminosity is only a factor of a few lower than that of the whole Galactic stellar halo. As a result of this large initial extent, the debris occupies a large volume in phase-space. The top panels show X-Y distribution at three different times. The black circle indicates the location of a "Solar neighborhood" sphere of 4 kpc radius. In the bottom panels grey dots show all satellite particles whereas the black dots represent those inside the sphere, and reveal the presence of multiple streams in the system. Credits: Adapted from Gómez & Helmi (2010), their Figs. 4 and 5.
∆ 6 w ∼ ∆ 3 x∆ 3 v and since ∆ 3 x (the spatial extent covered by the debris or 1/ρ) grows in time, this means that ∆ 3 v decreases with time locally as shown by see Buckley et al. (2019) for a slightly different and interesting application of these concepts). This implies that in a given location in the Galaxy, one may have many different moving groups sharing a common origin, as clearly apparent in Figure 4f which depicts a phase-space slice of stars in a simulation of a relatively massive accreted satellite.
From these considerations, it transpires that to detect each of the predicted moving groups large samples of stars are needed with accurate kinematics.  estimated analytically (this was later confirmed using cosmological simulations by Helmi et al. 2003) that if the whole stellar halo had been built via mergers, approximately 500 streams would be expected in the halo near the Sun (independently of whether 10 or 100 galaxies had been accreted). Given that the velocity dispersion of the halo is ∼ 100 km/s, the velocity resolution required would be 100/(500) 1/3 ∼ 13 km/s, and the sample size needed would have to contain at least as many as 5000 halo stars to yield on average 10 stars per stream. These estimates have nearly been met by Gaia DR2. Of course, higher precision and larger numbers of tracers would be necessary to go beyond simple detection of granularity (Gould 2003) to the full characterization of the streams and their parent objects.
As described above, debris originating in a single galaxy is thus expected to have similar integrals of motion (which includes of course, the adiabatic invariants). This has led to the search for ancient accretion events by looking for lumpiness in a space of "Integrals of Motion" (IoM). The first application of this method was by  which led to the discovery of the Helmi streams. Then a proof of concept of what would be possible with a mission like Gaia was given in Helmi & de Zeeuw (2000). Diagrams of E vs. Lz or Lz vs. L ⊥ (where L 2 ⊥ = L 2 − L 2 z , acts as a proxy for a third integral) are now being widely used to establish Galactic accretion history. The advantage of using IoM is that all the individual streams (or wraps) of a single object fold as it were into defining a single clump (compare for example, the ensemble of grey points in panels c and e to panel b in Figure 4). Therefore, the precision required on the measurements is less demanding and the signal of a clump in IoM is higher, because the number of stars in a clump is the total number of stars from each of the streams from a given object added together.
There are two caveats however. In an ideal case, energy or other integrals of motion would be conserved. However, the gravitational potential in which the streams have evolved must have changed with time, implying that this is not exactly true. Nonetheless, for example Gómez et al. (2013), Simpson et al. (2019) have shown that substructure is still present in these spaces, even in simulations of the full hierarchical assembly of the halo. Actions being adiabatic invariants, are better conserved although more difficult to compute (but see e.g. Sanders & Binney 2016). Thus far however, there has not been a real need to resort to them for the identification of merger debris. The likely reason is that the volumes probed so far by data with full phase-space information (6D) are sufficiently small that in the expression E = 1/2v 2 + Φ(x), the potential term is approximately constant, i.e. Φ(xsun) = Φ0, and so time variation, or even limited knowledge of the exact form of the potential has not been a limiting factor. The situation will change as we begin to explore beyond the solar neighborhood, especially with Gaia DR3 and subsequent Gaia data releases.
On the other hand, only if all the stars from a given accreted system would be mapped, the clump defined would be fully smooth (in the absence of dynamical friction). As discussed above, when we observe locally we typically only probe portions of debris streams. This implies that we expect substructure to be present within a clump associated to a given object in IoM space when using spatially localized samples of stars. This is clearly seen in Figure 4c, where the grey particles denote all the stars from the system (independent of their final location within the host) and those in black only those inside a small volume (indicated by the circle in the top right panel of Figure 4). Substructure in IoM may also appear if the system is very massive, and thus suffered dynamical friction. In that case, the orbit will have changed with time, and material lost early can be on significantly different orbits than that lost later.
Individual (portions of) streams are particularly apparent in frequency space, as can be seen in Figure 4a. This is because the individual streams each have their own characteristic frequency (which defines their phase along the orbit, see McMillan & Binney 2008, Gómez & Helmi 2010. The regular pattern seen in Figure 4a depends on the time of accretion of the system since ∆Θ = ∆Ωt, where ∆Θ represents the difference in phase, and ∆Ω would be the separation between neighboring clumps, i.e. a characteristic scale in frequency space. Therefore, since the stars plotted in this figure have all roughly the same location but differ in phase by ∆Θ ∼ 2πn (with n an integer), this implies that t could be inferred by applying a Fourier analysis, provided enough stars are found in each stream (Gómez & Helmi 2010). It turns out frequency space is also useful for constraining the mass growth or time variation of the gravitational potential as the characteristic regular pattern becomes distorted depending on how the system has evolved (Buist & Helmi 2015. It may be possible to measure these effects using samples of nearby main sequence halo stars, as these stars are numerous and their velocities and distance estimates may be more accurate because of their relative proximity.

Generalities
Mergers play a key role in the hierarchical cosmological paradigm. This is, after all, by and large the way that galaxies build up their dynamical mass, i.e. their dark halos (Wang et al. 2011). Therefore, tracking mergers becomes of great importance in the goal of unraveling the build up of galactic systems. The only way we have to track past mergers over long timescales is by resorting to stars. This is why the stellar halo of the Galaxy could be considered the component to disentangle the merger history of the Galaxy. It is here where disrupted galaxies, cannibalized by the Milky Way, will most likely have deposited their debris. Some debris may be deposited in the thick disk by satellites on low inclination orbits (Abadi et al. 2003). It is also a place where we may find heated stars from the disk, i.e. from those present at the time of the mergers and which were perturbed on to hotter orbits (Zolotov et al. 2009, Tissera et al. 2013. Most of the mass in the inner regions of Milky Way-like dark halos are predicted to originate in a few massive progenitors (Helmi et al. 2002, Wang et al. 2011, implying that these will have hosted sizable luminous galaxies (Cooper et al. 2010). Therefore most of the information regarding these mergers will be traceable in the stellar halo.
Not only is the stellar halo interesting from the point of view of the merger history, but as mentioned earlier also because it contains (some of) the oldest stars and the most metal-poor ones (possibly together with the bulge). This is not necessarily a coincidence. The existence of a mass-metallicity relation for galaxies, implies that the proto-Milky Way will generally have been the more massive object in its cosmic neighborhood. This implies that accreted galaxies will have been less massive than the proto-Milky Way and hence on average, more metal-poor than the disk. Since these objects deposit debris in the stellar halo, it is natural for it to have a lower metallicity. (Of course, this shifts the question to a different one, which is understanding why and how such a mass-metallicity relation arises, see e.g. Tremonti et al. 2004). Since there is also a correlation between mass and star formation rate (SFR hereafter), and even though the first stars to form in the Galaxy might have been very metal-poor (or Pop III), the ISM of the proto-galaxy likely was enriched more quickly because of its high SFR, reaching quite fast a higher overall metallicity, as observed for example in the Galactic bulge/bar region (see e.g. Matteucci et al. 2018).
Understanding the age distribution is trickier, no less because of the fewer precise constraints. However, there is a simple explanation for why the halo should generally be older than the thin disk. Since mergers were much more frequent in the past, it is only after the major epoch of merger activity that a thin disk could grow to its full current extent. The concordance model predicts that the first stars will form in the highest density peaks, which will collapse first and which are typically associated to the more massive objects at later times (e.g. Diemand et al. 2005). This would mean that the first stars in our cosmic environments ought to have formed in the proto-Milky Way. Cosmological simulations suggest these first stars are likely part of the bulge or inner spheroid (White & Springel 2000, Tumlinson 2010, Starkenburg et al. 2017, El-Badry et al. 2018, where the outer halo would be slightly younger. Thus a slight age gradient (remember we are still discussing the epoch before the thin disk as we know it was in place) could arise from the fact that lower mass objects typically form their first stars a bit later. Later accreted objects would also have continued forming stars longer, and so contributed to the trend (Carollo et al. 2018). An age gradient was what Searle & Zinn (1978, SZ, hereafter) discovered when studying the age distribution of halo globular clusters, and what led to the SZ-fragments model of the formation of the halo. Not only outer globular clusters are younger but this is also apparent in recent studies of blue horizontal branch stars (Santucci et al. 2015). Note however, that the age distribution as well as the metallicity, particularly of the outer halo could be rather patchy, and could depend on the specifics of the merger history (e.g orbits, time of infall) and mass spectrum of accreted objects (e.g. Font et al. 2006).
In summary, because the stellar halo contains in proportion more pristine stars, it gives us a window into the physical conditions present in the early universe (e.g. Frebel & Norris 2015), and also on the early phases of the assembly of the Milky Way. Hence its relevance in a cosmological context.

State of the art / Most recent discoveries
Our knowledge of the Galactic halo has increased greatly in the last 20 years. Relatively deep wide-field photometric surveys such as SDSS (York et al. 2000) and PanSTARRS 13 (Chambers et al. 2016), and DES 14 (Abbott et al. 2018) more recently, have revealed large overdensities on the sky and many narrow streams (Bernard et al. 2016, Shipp et al. 2018). These are direct testimony of accretion events that have built-up the (outer) halo, as discussed in e.g. the reviews by Belokurov (2013) and Grillmair & Carlin (2016), as well as other articles in the book edited by Newberg & Carlin (2016).
The second data release of the Gaia mission is, on the other hand, currently driving a true revolution in our understanding of the (inner) Galactic halo. This might have been expected because of the need for full phase-space coordinates for large samples of stars to pin down formation history discussed in Sec. 3.3. What was unexpected perhaps was the discovery that a large fraction of the halo near the Sun appears to be constituted by the debris from a single object, named Gaia-Enceladus ). This object is sometimes referred to as Gaia sausage because of its kinematic signature . The other very important contributor in the vicinity of the Sun to stars on halo-like orbits is the (tail of the) Milky Way thick disk (Gaia Collaboration et al. 2018a, as can be seen in the left panel of Figure 2. These (proto-)thick disk stars have likely been dynamically heated during the merger with Gaia-Enceladus , Di Matteo et al. 2019. We elaborate on these points below.

Gaia-Enceladus.
Although the presence of stars with metallicities typical of the thick disk but with halo-like kinematics had been reported before Gaia DR2 (most recently in e.g. Bonaca et al. 2017), the distinction in the kinematics had not been so clearly seen until DR2, as can be appreciated from the comparison between the left and right panels of Figure 5, and by inspection of Figure 6 compared to the left panel of Figure 2. For stars within 2.5 kpc from the Sun and with |V − VLSR| > 200 km/s, i.e. traditionally the regime of the halo, approximately 44% of the stars are in the "hot" thick disk region (200 < |V − VLSR| < 250 km/s), while a large fraction of those remaining (between 60% and 80% depending on the exact definition) are in the elongated structure that is due to Gaia-Enceladus and indicated in the right panel of Figure 5     Map of the kinematics of halo stars in the sample used by Belokurov et al. (2018). This was obtained from the cross-match of the positions in Gaia DR1 and SDSS and use of the long time baseline to derive proper motions. The top panels show the distribution for stars in different metallicity bins, while in the bottom panels the residuals resulting from a Gaussian mixture model are plotted, with the contribution of 2 and sometimes 3 subcomponents as indicated by the ellipses. Although a complex kinematic distribution can be retrieved statistically, comparison to the left panel of Figure 2 using Gaia DR2 data reveals how the striking increase in quality of this dataset leads to a true distinction of the various components, as was also noted in Figure 5. The V R asymmetry seen in the rightmost bottom row panel for the faster moving component is also apparent in Figure 2, and likely the result of the impact of the Galactic bar on the kinematics of these stars, as also reported by Antoja et al. (2015) for the canonical thick disk stars. Credits: Reproduced from Belokurov et al. (2018), top panels of their Figs. 2 and 3.
These findings link to what was arguably one of the first stunning surprises on the halo in Gaia DR2: namely the color-(absolute) magnitude diagram of stars with "halo"-like kinematics (i.e. selected to have tangential velocities VT 200 km/s) and presented in Gaia Collaboration et al. (2018a) revealed the presence of two clearly distinct sequences, as shown in the top panel of Figure 7. These well-defined sequences point to the presence of distinct stellar populations (i.e. with different ages and metallicities), and are evocative of a "dual" halo (see Carollo et al. 2007, and discussed in some detail in Sec. 4.2.2). Gaia Collaboration et al. (2018a) tentatively suggested that the older and more metal-poor sequence corresponded to low α-abundance stars on retrograde orbits first reported by Nissen & Schuster (2010, 2011. On the other hand, Koppelman et al. (2018) demonstrated that this sequence was dominated by the large kinematic structure (or "blob" as it was referred to by the authors), seen in the right panel of Figure 5.
Driven by these findings, and by the fact that the mean motion of the stars in the kinematic structure was slightly retrograde (as appreciated from the location of the red vertical line in the right panel of where there is overlap with thick disk (which by and large must have formed in-situ), this immediately implies that the stars formed in a different system than the thick disk, as [α/Fe] will generally decrease as [Fe/H] increases. This means that these stars must have been accreted. Furthermore because the majority of the stars in the nearby halo are part of the "blob" (if they are not in the tail of the thick disk), this implies that a large fraction of the halo near the Sun has been accreted. The accreted system is what has been called Gaia-Enceladus. The presence of a significant α-poor (low Mg) chemical sequence was first reported by Hayes et al. (2018), who used data from the APOGEE survey. It was however hard for these authors to determine to which component the stars in this sequence belonged

Gaia-Enceladus/Gaia-Sausage
The prominence of Gaia-Enceladus had in fact been noticed by Belokurov et al. (2018) using an unpublished catalogue of proper motions obtained by combining SDSS and Gaia DR1 astrometric positions (also used in Deason et al. 2017). Although the separation between the thick disk and the halo is less sharp because of the lower astrometric precision (as can be seen from Figure 6), the authors found (after performing a Gaussian mixture model), that a high fraction of the halo stars had very large radial motions. Through comparison to zoom-in cosmological simulations, this significant radial anisotropy was interpreted as implying that the halo stars originated in a significant merger the Galaxy experienced between redshift 1 and 3. In the spirit of what was known before Gaia DR2, this however was not the only possible interpretation, particularly because the lower quality of the proper motions did not reveal a retrograde mean signal (and hence somewhat "abnormal") in the multi-Gaussian component decomposition, and chemical abundance information (in particular the sequence of [α/Fe]) was not used in the study. As the authors themselves acknowledge in their paper, perhaps this blob or "sausage" structure as it was called (see Myeong et al. 2018b, and the available versions on the ArXiv), was the result of a monolithic-like collapse of the kind proposed by Eggen et al. (1962), in the traditional model of formation of an in-situ halo. The result of such a collapse would likely put the stars formed on radially biased orbits. Now with the knowledge provided by Gaia DR2 data in combination with that of the APOGEE survey, revealing respectively the retrograde mean motion of the halo and the distinct chemical sequence defined by the majority of its stars, there is absolutely no question that Belokurov et al. (2018) had seen Gaia-Enceladus' mark in the kinematics of halo stars, and interpreted with remarkable insight correctly its accreted origin.
to because of the lack of proper motion information (their study was carried out just a few months before Gaia DR2) or the magnitude of this population. Nonetheless from the line-of-sight velocities, Hayes et al. (2018) concluded that the stars had "halo"-like kinematics and because of the characteristics of the sequence that they likely represented an accreted population. These considerations led Fernández-Alvar et al. (2018) to fit a chemical evolution model. These authors showed that the sequence could be reproduced if the stars had formed in a system with an average SFR of 0.3 M /yr over a period of 2 Gyr or so. By integrating this SFR, Helmi et al. (2018) subsequently estimated a stellar mass of ∼ 6 × 10 8 M for Gaia-Enceladus. Since the existence of a sequence had been known for some years since Nissen & Schuster (2010), there were studies in literature that had compared the ages of the stars in the sequence to those in the thick disk sequence in the metallicity range −1 [Fe/H] −0.6 (Schuster et al. 2012, Hawkins et al. 2014, more recently Vincenzo et al. 2019. The low-α stars were found to be younger than those in the thick disk. Since these are the most metal-rich stars and likely formed in Gaia-Enceladus before it was fully disrupted, this dates also the time of the merger, and at the same time it demonstrates that a disk was already in place in the Milky Way at the time, roughly 10 Gyr ago (see Gallart et al. 2019, and Sec. 5).
The details of the merger have still to be pinned down. For example, a coarse comparison of the Gaia kinematical data to existing simulations of the merger of a disk galaxy with a massive disky satellite by Villalobos & Helmi (2008, 2009, as shown in Figure 8, suggests that the merger was counter-rotating because the mean rotational motion of associated halo stars in the solar vicinity is (slightly) retrograde. The presence of specific features in velocity space, namely the arc seen in the right panel of Figure 5, and their resemblance to those seen in the simulations, support the retrograde infalling direction and also suggest that the merger's inclination was initially approximately 30 o . However other configurations might also be possible, as similar characteristics are found for a merger that is coplanar but where the accreted object is spinning in the opposite sense than the host as shown by Bignone et al. (2019) using the EAGLE cosmological simulations suite.

Figure 8
Toomre diagrams obtained from a set of simulations aimed to study the formation of the Galactic thick disk via a 20% mass ratio merger (Villalobos & Helmi 2008). In these simulations, a merger between a disky satellite and a host disk is modeled for different orbital configurations. In the left panel the orbit is prograde, while on the right panel it is retrograde. In both cases the initial orbital inclination is 30 o . The star particles plotted are located inside a solar neighbourhood-like volume, 4 Gyr after infall, with black and blue corresponding respectively to those from the host and from the satellite. Some differences are apparent in the kinematic properties of the merger product, one of the most noticeable being the presence of an arc at high rotational velocities, positive for the prograde and negative for the retrograde case. This arc is composed by star particles lost early during the merger before the object had fully sunk in via dynamical friction (see R. Bos, BSc thesis, Univ. of Groningen, http://fse.studenttheses.ub.rug.nl/20160/). The arc on the right hand side panel is very reminiscent of that seen in the Gaia DR2 data and shown in the right panel of Figure 5. Also Bonaca et al. (2017, their Fig. 6) report the presence of a similar structure in the LATTE cosmological simulation of a Milky Way-like galaxy, but its origin is not discussed in the paper. Helmi et al. (2018) suggest a mass ratio of ∼ 4 : 1 for the merger (which has been confirmed by Gallart et al. 2019) on the basis of the estimated stellar mass of Gaia-Enceladus, the expected stellar mass-to-halo mass ratio for objects of this size, and assuming a 10 10 M thick disk' stellar mass present at the time. Other authors have suggested stellar mass of 10 9 M up to 5×10 9 M (Mackereth et al. 2019. Clearly these estimates are uncertain but all pin-point to a significant merger, with lower masses (as proposed by Donlon et al. 2019) being inconsistent with the chemical abundances of the stars (i.e. they would violate the mass-metallicity relation) 15 . The exact configuration of the merger will have to be pinned down with more detailed modeling and by probing the 3D motions and properties of stars well beyond the solar neighborhood.
In summary, the halo near the Sun (if not in a "hot" thick disk) is dominated by debris

Mean rotational velocity of Gaia-Enceladus debris
Determining the mean motion of the debris from Gaia-Enceladus is important as it reveals properties regarding the initial configuration of merger. Support for a retrograde encounter is apparent from the comparison of the simulations shown in Figure 8 to the data shown the right panel of Figure 5, as well as from the location of the vertical line (denoting null rotation) in the latter figure. To quantify this further, we select stars on the [α/Fe]-poor sequence, for [Fe/H] ≥ −1.3, i.e. the region dominated by Gaia-Enceladus debris (see Figure 7). By imposing a 20% relative parallax (or distance) errors quality cut, using distances determined either via 1/ or from the Bayesian method by McMillan (2018), a selection is made that satisfies −500 ≤ Lz ≤ 500 kpc km/s. The resulting sample is small but within 1 kpc from the Sun, v φ = −21.1 ± 1.8 km/s (6 stars) and within 2 kpc v φ = −16.1 ± 2.8 km/s (23 stars), where the uncertainty has been estimated from re-sampling the individual velocity errors. For the stars within 1 kpc, these errors are all smaller than 7.2 km/s (the mean is 3.3 km/s), while for those within 2 kpc, the mean error is 10.8 km/s. The signal is therefore weak but it is robust to measurement uncertainties (statistical and systematic). This is because within 1 kpc from the Sun, the parallax is 1 mas, the relative errors are small ( 20%) for the bright stars in the Gaia RVS sample (Gaia Collaboration et al. 2018d), and the systematic parallax bias , Arenou et al. 2018) is negligible within this parallax range (even if as large as 0.05 mas). A similar amplitude mean retrograde motion has been reported by Mackereth et al. (2019) using APOGEE distances.
from Gaia-Enceladus, a very massive object that was accreted 10 Gyr ago. As such it most likely represents the last significant merger that the Milky Way experienced. Probably after this was completed, the (current) Milky Way thin disk could start a more quiescent growth phase, and this would be consistent with the ages of its oldest stars. Other large mergers the Milky Way has experienced since then include that with the Sagittarius dwarf, but because of the late the time of infall (∼ 8 Gyr ago) and mass ratio of 1 -5% (M ∼ 5 × 10 8 M , see for example Dierickx & Loeb 2017, Fardal et al. 2019), this has resulted in a less dramatic impact. Even the ongoing merger with the Large Magellanic Cloud is less important (with a mass ratio is probably ∼10%), although in both cases we do see their effect on the disk of our Galaxy, in the form of phase-space spirals 16 (Antoja et al. 2018), and waves (vertical and radial oscillations, see e.g. Laporte et al. 2018a.
Evidence with hindsight. Besides Belokurov et al. (2018), several other authors had in fact come across traces of Gaia-Enceladus without knowing. We had a rather fragmented view of the inner halo until Gaia DR2 because it was based on small samples and hence our knowledge was "patchy" or "incomplete". For example, the sample with the low-α sequence of Nissen & Schuster (2010), contained less than 100 stars (which nonetheless was a significant increase in comparison to the original discovery paper, which had just 13 halo and 16 thick disk stars, see Nissen & Schuster 1997). Although these authors were able to put forward a scenario rather similar to what has now been revealed with Gaia DR2 and APOGEE, it is not until recently that the realization came that most of the halo follows the Nissen & Schuster low-α sequence, thanks to the work of Hayes et al. (2018) on the chemistry and Koppelman et al. (2018) on the kinematics (see also Haywood et al. 2018). Possibly the reason for this is that the halo near the Sun is known to peak at a metallicity of [Fe/H] ∼ −1.6 and have [α/Fe] ∼ 0.2 − 0.4 around this iron abundance. It is only at higher metallicity that the sequence by Gaia-Enceladus becomes clearly separate and has lower [α/Fe] than the thick disk as seen from the lower panel of Figure 7. Other hints of this structure were present in the  study based on Hipparcos data. These authors identified a concentration of stars for [Fe/H] ∼ −1.7 with very high eccentricities (e ∼ 0.9), but they interpreted it as being stars formed during an ELS-like collapse in the early Milky Way. Brook et al. (2003) compared the Beers et al. (2000) sample to a cosmological hydrodynamical simulation of the formation of a Milky Way-like galaxy. In this comparison, the authors identified a clump of stars with retrograde motions which they proposed could be satellite debris. The presence of clumps of stars with retrograde motion goes back even further, to the work of e.g. Majewski et al. (1996) and Carney et al. (1996) and more recently Kinman et al. (2007). The globular cluster Omega Cen became the natural culprit, because it was known to have a retrograde orbit and because of its peculiar multiple populations, which led to the suggestion that it could have been the nucleus of a disrupted dwarf galaxy (Dinescu 2002, Bekki & Freeman 2003, Meza et al. 2005 .

The dual halo. The two sequences revealed in the CMD of Gaia Collaboration et al.
(2018a) are a direct reminder of the suggestion of a "dual halo" in the Milky Way. As just discussed it has been possible to attribute the first, more metal-poor sequence largely to Gaia-Enceladus. On the other hand, the second sequence is populated by stars which kinematically are clearly part of the tail of the thick disk seen from the right panel of The idea of a dual halo was discussed already quite thoroughly by Norris (1994) (see also the nice historical introduction in Carollo et al. 2010). Norris (1994), like several other authors around the same time, found that at lower metallicities the retrograde component becomes more and more prominent. In his interpretation of the data, this dual halo would consist of an accreted component (a la Searle & Zinn 1978) and a contracted halo (a la ELS, Eggen et al. 1962).
We now know from the Gaia DR2 data that what Norris (1994) called "the contracted halo" is actually largely heated (proto)-thick disk. So it was indeed formed mostly in-situ, however the stars did not form during a collapse but in a disk that was heavily dynamically perturbed during a merger, as in the simulations of Zolotov et al. (2009), Tissera et al. (2013). On the other hand, the accreted component is predominantly debris from Gaia-Enceladus. This explains why the local inner halo was for a long time considered to be slightly prograde (since it is a mix of prograde thick disk and slightly retrograde Gaia-Enceladus, and was concealed as a single component because of the large velocity errors as shown in the left panel of Figure 5), and more metal-rich locally than in the outskirts. In comparison, the debris from Gaia-Enceladus dominates farther out (at several kpc away from the Sun and above the Galactic plane), implying that the outer halo consists of more metal-poor stars on retrograde orbits (at least in comparison to the proto-thick disk). This is also what Morrison et al. (2009) had found in a carefully constructed sample of nearby metal-poor red giant branch stars with accurate distances, and what the more recent study by Bonaca et al. (2017) using Gaia DR1 data (Gaia Collaboration et al. 2016) also revealed. Carollo et al. (2007Carollo et al. ( , 2010 and Beers et al. (2012, and references therein) based on SDSS and SEGUE data put the idea of a "dual halo" on firmer ground, despite concerns that the retrograde mean motion signal was driven by biases in the distances as suggested by Schönrich et al. (2011). The Gaia DR2 data has put the debate fully to rest since the retrograde signal is very strong and is seen also in the nearby halo, where the distances have been measured very precisely and accurately with trigonometric parallaxes (or with Bayesian methods, McMillan 2018). There is however, still some tension with the work of Carollo et al. (2007), because what they called outer halo peaked at a metallicity of [Fe/H] ∼ −2.2, lower than what is typical for Gaia-Enceladus stars. A possible way out would be to consider that large galaxies like Gaia-Enceladus probably had a metallicity gradient. As the galaxy was destroyed, some of the material in the outer regions would be put on more unbound orbits. This material would have lower metallicity and form part of the outer portions of the inner halo (i.e. we are referring here to inner halo as that within 20 kpc from the Sun). Another possibility could be an issue with the metallicity scale calibration, or perhaps a small bias in the sample selection, which because it is based mostly on the colors for main sequence turn-off stars, might preferentially select more metal-poor stars ).

More substructures in the inner halo.
4.2.3.1. Helmi streams. One of the first genuine accreted inner halo kinematic substructures discovered near the Sun were the Helmi streams 17 . Before their discovery, there had been reports of substructure (e.g. Majewski et al. 1994) and in particular, Eggen argued for the presence of streams in the halo but their reality was often questioned because of the methods used (because in some cases, for example, the very uncertain distances were enforced to be such as to match those expected for a kinematic group 18 ). The Helmi streams were identified in a dataset largely compiled by Chiba & Yoshii (1998), who used the Hipparcos dataset supplemented by line of sight velocities and distances from the literature. The streams are also readily apparent in Gaia DR2 (Gaia Collaboration et al. 2018c, Koppelman et al. 2018). A characteristic of the streams is that they have rather high z-velocities, which makes them easy to identify as they are well separated from the rest of the stars with the same rotational velocity (v φ ∼ 150 km/s). One of the streams has positive vz while the second (and better populated) has vz < 0. Koppelman et al. (2019b) have shown that the streams likely stem from a progenitor of M * ∼ 10 8 M . Seven globular clusters have been associated to this object, including e.g. NGC 5024 and NGC 5053 (Koppelman et al. 2019b, Massari et al. 2019. These clusters follow a welldefined age-metallicity relation and their number is consistent with that expected from the globular clusters' specific frequency relation given the progenitor's mass (see for example Zaritsky et al. 2016). The time of infall estimated by Koppelman et al. (2019b) is between 5 and 8 Gyr ago (see also Kepley et al. 2007) and it is largely driven by the asymmetry 17 This naming has been used informally in workshops and conferences for years and it appears in written print in Klement (2010).
18 It would nonetheless be an interesting exercise to revisit in a systematic fashion the Eggen moving groups with the data newly available, particularly from Gaia DR2. For example, Navarrete et al. (2015) have analysed the Kapteyn group and the chemical abundances of its stars and shown that they do not really constitute a physical unity. seen in the number of stars with positive and negative vz. It is somewhat puzzling that an object that is not so massive would have made it to the inner halo so recently. Since its orbit seems to lie close to a resonance (J. Hagen, PhD thesis, Univ. of Groningen, 2020), the inferred timescale could be underestimated as streams will then have spread out more slowly. Alternatively the progenitor system may have fallen in together with a heavier galaxy as part of a group (as the Magellanic Clouds, e.g. Sales et al. 2017).  There are now enough members of the Helmi streams known that is is possible to produce a well-populated CMD, which is very similar to those constructed for many decades for the dwarf spheroidal satellites of the Milky Way (see e.g. Tolstoy et al. 2009) as shown in the left panel of Figure 9. The information contained in this CMD should allow deriving a star formation history. In combination with the metallicity distribution function shown in the right panel of Figure 9, this should enable the reconstruction of the full chemical history of the system (in a similar manner as done for several dSph in e.g. Homma et al. 2015). Detailed elemental chemical abundances (such as those already measured by Roederer et al. 2010), might also help in this respect (see for example Lanfranchi & Matteucci 2010). It is intriguing that such an exercise of reconstructing the star formation and chemical enrichment history is now possible using nearby stars and for objects long-gone. This truly lies at the heart of what it means to carry out Galactic archaeology.

Other reports of substructures. Several other substructures have been reported
in the literature prior to Gaia DR2. Typically though, these structures were less conspicuous as they were based on Hipparcos/Tycho data in combination with ground-based radial velocities (see the reviews by e.g Klement 2010, Smith 2016. Kepley et al. (2007) and Morrison et al. (2009) report lumpiness in a sample of red giant branch stars with accurate distances, some of which can be attributed to Gaia-Enceladus and to the hot thick disk, but also "outliers" on very prograde orbits. More recently other authors used LAMOST and TGAS (Liang et al. 2017 Koppelman et al. (2018) have been suggested to be part of a larger substructure associated to the debris of a dwarf galaxy dubbed as "Sequoia" . The associated stars have more retrograde motion and a lower metallicity than Gaia-Enceladus as can be seen from Figure 10, and they are on a less eccentric orbit. Myeong et al. (2019) suggest that Omega Cen could have been Sequoia's nuclear star cluster, and that several other clusters, including the recently discovered FSR 1758 (Barbá et al. 2019) would also be associated. On the other hand, Massari et al. (2019) argues that it is more likely that Omega Cen is associated with Gaia-Enceladus, given its location in the age-metallicity diagram of the Galactic globular clusters. As can be seen the right panel of Figure 10, the extent of Sequoia in the space of energy and Lz if the two clusters would be associated to it, is as large as that of Gaia-Enceladus, implying that it would need to have been as massive (see Koppelman et al. 2019a, for details). This however would be in tension with its lower metallicity.

Sequoia and Thamnos. Two of the clumps identified in
The presence of additional (besides Gaia-Enceladus) debris with retrograde motions put forward by Myeong  On the other hand, Matsuno et al. (2019), using chemical abundances from the SAGA database 19 (see Suda et al. 2008, and subsequent papers), showed that the high-energy more retrograde stars have lower [Mg/Fe] at [Fe/H] −1.6 than Gaia-Enceladus stars and especially, different [Na/Fe]. The authors further argue that the other elemental abundances are too different from those measured for Omega Cen for the cluster (which has a higher binding energy) to be related.
Finally, Koppelman et al. (2019a) report evidence that the lower energy region occupied by stars tentatively associated to Sequoia, is likely part of a different separate structure, which the authors name "Thamnos" and indicated in Figure 10. Thamnos stars define the separate clump to the lower left of the region occupied by Gaia-Enceladus (vy ∼ −150 km/s  and v ⊥ < 150 km/s). Their progenitor would be a small galaxy accreted on a low inclination orbit, and this is consistent with the lower mean metallicity of its stars. Thamnos, which in ancient greek means "shrubs", overlaps also in part with substructures that have been reported earlier in Helmi et al. (2017), Koppelman et al. (2018). The majority of the low eccentricity stars with high [Mg/Fe], low metallicity and retrograde motions identified by Mackereth et al. (2019) overlap with those from "Thamnos".
It should also be noted that Figure 10 reveals that some of the stars from Sequoia overlap directly with the arc-like structure assigned to Gaia-Enceladus on the basis of the resemblance to numerical simulations as discussed before. If Gaia-Enceladus had a metallicity gradient, which might be expected given it was of comparable size to the Large Magellanic Cloud, then one might expect the outskirts to be more metal poor and part of the material lost first, to be on a more a retrograde orbit, as found for some Sequoia stars.
The above discussion serves to demonstrate the high-state of flux of the field and that it is currently in some turmoil. It evidences that it is very difficult, given the available data and methods, to pin down the origin of the different clumps and overdensities being discovered in the halo near the Sun. This is particularly manifested in the analysis of Yuan et al. (2019) who identified more than 50 small clumps using neural network based clustering on the very metal-poor sample of stars from LAMOST. Nonetheless, it is possible that very detailed modeling and chemical analysis of a large sample of stars will help shed light on their true nature. At the moment, the modeling is sketchy, carried out in idealized (timeindependent or non-cosmological) configurations, often hydrodynamics and star formation (as well as chemical enrichment) are not followed, and the numerical resolution (especially

The very retrograde halo
Put together, a possible scenario to explain the findings reported in this Section would be the existence of (a bonsai) Sequoia as embodied by a clump of high-energy retrograde stars that is truly distinct from Gaia-Enceladus (because of the differing Na/Fe abundances reported in Matsuno et al. 2019), but whose extent is smaller than that originally proposed by Myeong et al. (2019). It is not inconceivable that the region of the arc-like structure contains a mix of debris from two different objects, stars from the outskirts of Gaia-Enceladus and from the Sequoia galaxy. An alternative possibility would be that we are seeing only stars from the outskirts of Gaia-Enceladus and that these followed a somewhat different chemical enrichment path. On the other hand, the lower energy region, with stars of low eccentricity and high [Mg/Fe] identified by Mackereth et al. (2019) and Koppelman et al. (2019a), would be associated to "Thamnos".
in terms of particle number) is typically too low to do a proper comparison to the features seen in the data. Furthermore, the samples with detailed and precise chemical abundances of halo stars, even for those that are bright and nearby, are relatively small. Finally, more robust statistical, or probabilistic tools are needed to assign stars to overdensities, as these very likely have some degree of overlap in phase-space.

Substructures in the outer halo.
At larger distances from the Sun, i.e. in the outer halo, many overdensities have been uncovered, especially with photometric wide-field surveys such as SDSS. Examples are the Sagittarius streams (Ivezić et al. 2000, Hercules Aquila Cloud (Belokurov et al. 2007) and the Virgo overdensity (Vivas et al. 2004, Jurić et al. 2008. At this point the nature of these high latitude structures (beyond the Sagittarius streams) and the link to those in the solar vicinity discussed in previous sections remains unknown, although Simion et al. (2019) have proposed that both Hercules-Aquila and Virgo are related to Gaia-Enceladus. To pin this down requires large samples of stars with accurate full phase-space information, some of which will become available in future Gaia data releases combined with data from large spectroscopic surveys such as WEAVE 20 (Dalton et al. 2016 Other known overdensities located closer to the Galactic plane, are the Monoceros ring (Newberg et al. 2002), and the feathers at slightly higher latitudes (Grillmair 2011, Bernard et al. 2016, seen in SDSS andin PanSTARRS, respectively). These low latitude structures are likely the result of the response of the disk to a massive perturber, such as Sagittarius, and perhaps also contain some satellite debris (Sheffield et al. 2018, Laporte et al. 2018b. Besides large overdensities, many thin streams crisscrossing the halo at larger distances have been uncovered (see Grillmair & Carlin 2016, Shipp et al. 2018). These include for example, GD-1, Atlas, Orphan, Pal 5 and many more. Mateu et al. (2018) has also identified 14 candidate streams using RR Lyrae stars and produced a very interesting and useful interface, the GalStreams Footprint Library and Toolkit for Python 23 , to visualize and keep track of all known structures. An example of the distribution on the sky of currently known streams and spatial substructures is shown in Figure 11. Thin streams tell us the story of the destruction of less massive objects, i.e. smaller dwarf galaxies and globular clusters, and are particularly useful for constraining the Galactic potential and the distribution of mass at large radii. An exciting development is the detection of substructure, gaps and overdensities in these streams, potentially revealing the presence of (dark) satellites orbiting the halo of the Milky Way .  Sky distribution of currently known spatially coherent streams and overdensities (indicated as regions delimited by dashed lines, and in boldface in the inset) produced using the GalStreams package by Mateu et al. (2018). Credits: C. Mateu and E. Balbinot.
The Streamfinder algorithm developed by  has also allowed the discovery of less distant thin streams. These include Gaia-1 and Gaia-2 , and the stream from Omega Cen itself (Ibata et al. 2019a). Streamfinder works by randomly sampling radial velocities (which have not been measured for the majority of the stars in Gaia DR2), while making use of the proper motion and photometric information of stars (i.e. for a given color there are at most 3 possible absolute magnitudes, see for details Ibata et al. 2019b). From these tentative phase-space coordinates, orbits are integrated in a Galactic potential and if stars can be found along the orbit on a "spaghetti" or stream-like structure, then a stream is identified (after proper statistical assessment and comparison to a suitable background). Watkins et al. (2009), Deason et al. (2013 have reported the presence of a break in the density profile of the Galactic stellar halo, at a distance of ∼ 20 − 25 kpc from the Galactic center. Deason et al. (2018) demonstrated that is likely marking the orbital turning-points (i.e. a shell), of the debris from Gaia-Enceladus'. These authors integrated the orbits of the stars belonging to the Sausage (Gaia-Enceladus), and showed that their apocenters lie at roughly this distance. Several shells might be expected from the debris, possibly at different distances because the object must have experienced significant dynamical friction. This break manifests itself also in a change in the shape of the velocity ellipsoid of halo stars around this distance, as reported in Kafle et al. (2014).

Link between inner and outer halo (sub)structures. Several authors including
In an attempt to make the link between the outer and the inner stellar halo even stronger, we have followed a similar approach to that of Deason et al. (2018), and have integrated the orbits of halo stars located within 1 kpc from the Sun. Since we expect the inner halo to be well mixed, these stars' trajectories should give us a broader view of the 3D (spatial) structure of the halo for example. This is much in the same way as orbits can be seen as the building blocks of a galaxy in Schwarzschild's modeling (Schwarzschild 1979) and can be used to reproduce for example their light profile.
The top panel of Figure 12 shows the orbits of all stars from Gaia DR2, located within 1 kpc from the Sun, with halo-like kinematics, defined as having |V − VLSR| > 210 km/s. In this case, we have used the distances from McMillan (2018). We have used their current positions and velocities (corrected for the Solar motion and the Local Standard of Rest) as initial conditions for the orbital integrations. These were performed in the MilkyWay Gala potential (Price-Whelan 2017, Price-Whelan et al. 2017), which contains an NFW halo, a nucleus, and a disk and bulge (Bovy 2015). The integration covers 2 Gyr. We have plotted each orbit sampled in 0.1 Gyr intervals, color-coded by their distance from the Sun. The top panel of Figure 12 reveals the presence of many overdensities, boxy shapes and some sharp edges, some of which bear some resemblance to those seen in the distribution of RR Lyrae in PanSTARRS by Sesar et al. (2017), and shown in the top panel of Figure 13.
The bottom four panels of Figure 12 show the orbits of stars now separated by the progenitor they are presumably associated with according to Koppelman et al. (2019a) roughly following what is shown in Figure 10). Comparison of the panels reveals that these four objects do indeed have different orbital properties as they have deposited their debris in different spatial configurations. For example, the debris from Gaia-Enceladus follows a symmetric configuration with respect to b = 0 o , with sharp boundaries at l ∼ ±120 o . The limited extent in longitude is similar to that reported in Helmi et al. (2018) and shown in the middle panel of Figure 13. This map shows the sky distribution of stars selected to be part of Gaia-Enceladus (with relative parallax errors < 20% and with a rather generous cut on Lz towards the retrograde halo, and a sharper more conservative one towards prograde moving stars to avoid contamination from the hot thick disk). The remarkable difference with the sky map obtained via the orbit integration is the degree of incompleteness and the effect of selection biases particularly for distant stars (with 0.2 mas) that appear as a result of the imposed parallax quality cut. Incompleteness also affects the RR Lyrae map of  shown in the bottom panel of the same figure, although in this case it is due to poorer time-sampling and hence characterization of their light curves Figure 12 Sky projections of the orbits of stars with halo-like kinematics currently located within 1 kpc from the Sun (and with |Lz| ≥ 10 kpc km/s). The orbits have been integrated in a Milky Way-like potential for 2 Gyr, and each point in the diagram corresponds to a point on a star's orbit color-coded by its distance. In the top panel we show all the stars together, while in the four lower subpanels we have separated the stars in the structures identified thus far following the assignment of Koppelman et al. (2019a), which is similar to that schematically shown in Figure 10. The decrease in the density of points towards the Galactic center is not physical and is due to removal of stars ( 1.3%) with very radial orbits of low |Lz| angular momenta . Top panel: Distribution of RR Lyrae stars from PanSTARRS. Middle panel: Distribution of potential members of Gaia-Enceladus from Gaia DR2 (with −1500 < Lz < 150 km/s kpc) and with parallaxes > 0.1 mas. The color coding indicates parallax, with more distant stars in yellow. The white contour encompasses 90% of all stars in the 6D Gaia DR2 sample with 0.1 < < 0.2 mas and 20% relative parallax error. Tentative members of Gaia-Enceladus follow a similar distribution on the sky as all the distant stars implying that the parallax quality cut introduces a selection bias: distant stars outside of these contours could be missed because these regions of the sky have been visited less frequently by Gaia. in Gaia DR2. This explains the significant asymmetries with respect to the Galactic plane seen in middle and bottom panels of Figure 13, which would not be consistent with the estimated time of accretion and the size of the progenitor, as its debris is expected to have fully phase-mixed. This is indeed what the velocities of nearby stars from Gaia-Enceladus predict, as shown in Figure 12 (as pointed out by H-W. Rix, private communication). On the other hand the distribution in longitude revealed in the studies of Helmi et al. (2018),  and also in Sesar et al. (2017) are not too dissimilar to that predicted for Gaia-Enceladus as seen in Figure 12. They all reveal a very centrally concentrated distribution of the debris. Figure 12 could also serve as guidance in searches for associated debris particularly at large radii. Furthermore, since the exact location of the features is sensitive to the gravitational potential, a comparison between the outcome of the orbit integrations and observational data (position in the sky, distance and kinematics of the various features) could be used to constrain better the distribution of mass in the Milky Way. Some of the features may even link to the overdensities discussed in the literature and already mapped at larger distances from the Sun, as a coarse comparison to Figure 11 suggests.

The thick/early disk
Not so many review articles have been written on the Galactic thick disk (but a good starting point is the introduction of Robin et al. 2014). This is likely because its reality (independent of that of the thin disk) has been highly debated throughout the years (Gilmore & Reid 1983, Bahcall & Soneira 1984 and also recently (Fuhrmann 2011, Bovy et al. 2012. Another likely reason is that sometimes conflicting answers regarding its properties have been obtained depending on the type of observational tool used to characterize its properties (e.g abundances, kinematics, star counts, see for example, de Jong et al. 2010, Cheng et al. 2012. This is discussed in Kawata & Chiappini (2016) and an insightful explanation is given in Minchev et al. (2015). We will not attempt to provide here a review but we will mostly focus on some observational facts and on recent discoveries, particularly in relation to the stellar halo, which help us understand at least in part the formation or evolutionary history of this component.

Overview of its properties
The thick disk was discovered through star counts by Gilmore & Reid (1983). These authors found an excess of stars at large heights above the plane, beyond what would be expected from a single exponential fit corresponding to the thin disk. The excess could be fit by invoking a second component following also an exponential functional form, but with a larger scale height. Subsequent work revealed that the stars in the thick disk had different kinematics, which although mostly rotationally supported, have lower rotational speeds (by about 30-50 km/s) and the velocity dispersions were higher than the thin disk. The first spectroscopic studies showed the thick disk to be more metal-poor than the thin disk and to be composed of stars which were older, see for example Sec. 4 of the extensive review by Gilmore et al. (1989).
More detailed high-resolution chemical elemental abundance studies demonstrated that thick disk stars organize themselves in a segregated sequence from that of the thin disk stars in the solar neighborhood, for example in [α/Fe] vs [Fe/H] (Gilmore et al. 1995). Several authors have recently provided definitive evidence that the sequences are truly separate, and hence that the two components are really physically distinct, as they are made up of stars that do not overlap in their properties (e.g. Adibekyan et al. 2011, Recio-Blanco et al. 2014, Hayden et al. 2015. Haywood et al. (2013) also showed that the stars in the thin and thick disks follow very tight and well defined tracks in [α/Fe] and [Fe/H] with age, with a break occurring at ∼ 8 − 9 Gyr, which marks the oldest stars present in the thin disk. These distributions display small scatter, a result that although based on a local sample can be extended beyond the solar vicinity since the orbits of the stars probe a relatively large radial range (from 2 -10 kpc from the Galactic center). This small scatter (which implies no radial gradient) can be explained if the (majority of) thick disk stars formed rather quickly in a massive gaseous disk, possibly supported by turbulence (Snaith et al. 2014).
The thick disk metallicity near the Sun peaks at [Fe/H] ∼ −0.5, and extends on the metal-rich side up to solar metallicity. It also has a very significant tail, which is often referred to as the metal-weak thick disk (Norris et al. 1985, Morrison et al. 1990, Morrison 1993, Beers et al. 2002. The stars associated to this tail are clearly visible as the blue/green data points with high [Mg/Fe] and [Fe/H] −1 in Figure 14, which is based on APOGEE data (Mackereth et al. 2019). This metal-weak thick disk could potentially be related to the very first disk or the oldest disk that was ever formed in the proto-Milky Way.

Formation paths
Typically four different scenarios are discussed in the literature for the formation of the thick disk (Gilmore et al. 1989, Robin et al. 2014). The traditional/oldest is via a minor merger onto a pre-existing disk, which leads to dynamical heating and the formation of a hotter but still rotation-supported component (Quinn et al. 1993). The accretion scenario is based on cosmological simulations which showed that if satellites are preferentially accreted from specific directions, this can lead to their debris being deposited in a planar configuration (Abadi et al. 2003). On this preferred plane, gas would later cool down and form the thin disk. The gas-rich scenario is inspired in cosmological hydrodynamical simulations that show that disks were highly turbulent and hotter in the past, partly because they were more gas rich and also because of the ongoing merger activity that prevented full settling (Brook et al. 2004, Bird et al. 2013. This is also what observations of high redshift disks appear to suggest (e.g. Bournaud et al. 2007). A last scenario is that of migration, that is stars from the inner (thin) disk have migrated with time to the outer regions of the disk.
Because of inside-out formation and metallicity gradients, these stars would be older and have different chemical composition. Schönrich & Binney (2009) were the first advocates of this model that have quantitatively explored its feasibility. Sales et al. (2009) proposed that one way to determine the dominant formation mechanism of the thick disk would be from the eccentricity distribution of its stars. These authors showed that the different paths discussed above led to different distributions, with radial migration changing only slightly the low eccentricities of the stars. On the other hand, a dry large minor merger would leave behind a distribution of stars with intermediate eccentricity (the heated disk) and a high eccentricity bump formed mostly by accreted stars. A comparison to data from RAVE and SDSS carried out later by Wilson et al. (2011) and by Dierickx et al. (2010) showed that the most likely path was through gas rich mergers, i.e. turbulent disks in which stars were forming during mergers (Brook et al. 2004). This interpretation and idea has been largely confirmed by the latest analyses based on Gaia DR2. As Gallart et al. (2019) have shown, the majority of the thick disk stars likely formed after/during the merger with Gaia-Enceladus, and not before. However, some fraction did form before, as in the dry manager scenario, although the predicted bump with high eccentricities associated to the accreted stars (Sales et al. 2009, Di Matteo et al. 2011 is not seen in the thick disk. In fact, these stars exist but now we know they make up a large fraction of the Galactic halo (i.e. this is Gaia-Enceladus debris). It is interesting that the connection between thick disk and halo had not been fully made until recently (although see Purcell et al. 2010, who using numerical simulations discussed this possibility).
Although it is probable that radial migration has played some role in the evolution of the thick disk, and that some fraction of the stars in the thick disk have an (inner) thin disk origin (see e.g. Adibekyan et al. 2011), it is now likely that the efficiency of this process was initially overestimated (as argued by Minchev et al. 2012). The evidence discussed above, and particularly the work of Gallart et al. (2019) supports a scenario in which a gas-rich disk experienced a merger with Gaia-Enceladus, where the stars already present were dynamically heated, and star formation was triggered (possibly in a starburst) leading to the formation of the bulk of the stars in the thick disk. This interpretation is based on what is shown in Figure 15. This figure reveals that the distribution of ages of stars in the thick disk "proper" (in black) peaks at ∼ 10 Gyr, while the stars in the "hot" thick disk (the red sequence in the top panel of Figure 7  This scenario is also consistent with the rather uniform distribution of [α/Fe] with age and radius discussed in Haywood et al. (2015), as such a merger likely triggered a global response of the whole disk. Although there was probably a large amount of radial migration during the merger as the disk became dynamically hotter, this migration (only "blurring" no "churning") would not have been due to internal mechanisms as proposed in Schönrich & Binney (2009), but externally induced.

Figure 15
Age distribution for stars with |Z| > 1 kpc and V T < 200 km/s (i.e. the thick disk, in black) and for stars on the red and blue "halo" sequences (with b > 30 deg and V T > 200 km/s) as derived from the analysis of Gaia DR2 photometric data. The stars on the blue and red sequences, which correspond to Gaia-Enceladus and to the "hot" thick disk respectively, are both old but have different colour because of their different metallicity. The thick disk "proper" (in black) is younger and more metal-rich. This figure reveals agreement with earlier work showing that the youngest stars in Gaia-Enceladus are 9-10 Gyr old (the tail for younger ages is likely contamination), and also shows that a starburst in the thick disk appears to have been triggered 10 Gyr ago. Probably this is the time of the closest encounter between the two interacting systems, after which Gaia-Enceladus was fully engulfed by the Milky Way.

Further insights on the early disk from chemistry and dynamics
The evidence accumulated so far and discussed above suggests that we may have identified the reality of the proto-disk as the metal-weak (or "hot") thick disk, and that this was present more than 10 Gyr ago. It is interesting that relatively low eccentricity stars have survived as such the dynamical impact of a massive minor merger 24 , allowing this disk to be traced back to lower and lower metallicities as Figure 14 shows. In fact in recent work, Sestito et al. (2019) on ultra metal-poor stars (with [Fe/H]< −4) have shown that a fair fraction of these (26%) actually are rotationally supported and have thick disk-like kinematics. These stars would potentially be tracing the most pristine disk in the Milky Way. It will be interesting to bridge the gap in metallicity between the metal-weak thick disk and the regime probed by the most metal-poor stars to trace the history of the very first disk-like component in our Milky Way. This is particularly relevant in the context of linking the Milky Way to studies of high-redshift disks (e.g. van Dokkum et al. 2013. Just like we have done for the halo, we can now put previous work on the thick disk in the recently gained context. For example, Gilmore et al. (2002) discovered an excess of stars towards the rotation fields (i.e. l ∼ 90 o , 270 o ), with lags of approximately 100 km/s (as well as a minor contribution from a retrograde component). They argue this is due to "shear" in the kinematics of the thick disk such that at higher latitudes (b ∼ 33 o , 45 o ), thick disk stars rotate more slowly than near the plane. They interpreted this as evidence for a shredded satellite, but instead (given the evidence we have discussed so far) it is likely they were seeing "kicked-out" thick disk stars, and a bit of debris from Gaia-Enceladus. They were nonetheless "Deciphering the Last Major Invasion of the Milky Way" as the title of their paper suggested. Further analysis of Wyse et al. (2006) and Kordopatis et al. (2013) towards other lines of sight confirmed their first results. Other evidence hinting at the dynamical consequences of a significant merger on the early disk, are the overdensities discovered by Larsen & Humphreys (1996), Larsen et al. (2011) suggesting that the thick disk may be triaxial. Such a configuration is not an uncommon end-product of simulations of disks experiencing a massive minor merger (Villalobos & Helmi 2008). In these simulations this shape is delineated only by some of the stars already present at the time of the merger, which as Gallart et al. (2019) argue, do not comprise the majority of the thick disk of the Milky Way.
Other evidence of "substructure" in the thick disk was put forward by Schuster et al. (2006), where two different groups of stars in the thick disk with different mean metallicities and mean rotational velocity were identified. This is in fact one of the key papers preceding the Nissen & Schuster (2010) discovery of the two sequences, since what Schuster et al. (2006) were seeing was in fact stars from Gaia-Enceladus and from the thick disk. Analysis of the Geneva-Copenhagen survey (Nordström et al. 2004) led Helmi et al. (2006 to also propose the presence of substructure in the region kinematically dominated by thick disk stars. What these authors demonstrate in a follow up paper ) is that there is a transition in the dynamical properties of stars at a metallicity of [Fe/H] ∼ −0.4. Below this value stars have a large range of eccentricities, while above it stars are only on low eccentricity orbits. There is also significantly more scatter in [α/Fe] below this [Fe/H] value, as if there were a mix of populations (as in fact shown by Mackereth et al. 2019, and reproduced here in Figure 14). Bearing differences in metallicity scales, this is the [Fe/H] at which the thick disk separates into stars formed before/after the merger with Gaia-Enceladus (Gallart et al. 2019). Thick disk stars below this value (i.e. more metal-poor) have been kicked-out on to more extreme orbits (and there may even be some contamination from Gaia-Enceladus, see Di Matteo et al. 2019). On the other hand, stars formed at higher metallicities formed and have stayed on (proper) thick disk-like orbits. This was not the original interpretation given by Helmi et al. (2014) (and follow up papers such as Stonkutė et al. 2012Stonkutė et al. , 2013Stonkutė et al. ,Ženovienė et al. 2015 where the features were attributed to the presence of merger debris, whereas we now believe the latter is a minor contributor and what is seen is simply the imprint of an important transition in the history of the disk (similar conclusions using different orbital parameters were in fact reached earlier by Liu & van de Ven 2012).
This is an important point as there has been some propensity to attribute substructure or overdensities to accretion events. As vehemently argued by Jean-Baptiste et al. (2017) this is not necessarily the case. A merger can induce asymmetries and substructure also in the populations formed in-situ (as shown already in e.g. Gómez & Helmi 2010, Gómez et al. 2012, for the thick disk simulations of Villalobos & Helmi 2008). On the other hand, asymmetries and substructures can also arise from internal dynamical processes by such as resonances with, for example, the Galactic bar and which is responsible for the Hercules stream (Dehnen 2000). Also stars in the thick disk are affected by the bar, as shown in Antoja et al. (2015), and as evidenced for example in the left panel of Figure 2, where the vR velocity distribution of the (hot) thick disk is clearly asymmetric in the same way as the thin disk, whose asymmetry is explained as being due to the bar.
The above discussion serves to stress that care is required in the interpretation of substructure. Nonetheless, if substructures can be proven to be related to mergers, this is interesting from the archaeological point of view. As just discussed such substructures can reveal both the response of the in-situ system (and hence contain information about its properties and the nature of the encounter), as well as on the accreted population.

Discussion
Gaia DR2 data supplemented with that from existing large spectroscopic surveys, have made it relatively straightforward to identify what was plausibly a very important milestone in Galactic history: the merger with Gaia-Enceladus. Stars from this galaxy distinguish themselves kinematically, having slightly retrograde and very eccentric orbital motions, and more patently via their chemical abundances. In particular, these stars define a tight chemical sequence in [α/Fe] vs [Fe/H], which is especially apparent for [Fe/H] −1.3. The sequence merges with that of the (metal-weak) thick disk for lower metallicities, and the values of [α/Fe] here depict significantly more scatter (see Figure 14). It seems improbable that the scatter is due to the overlap of stars from just these two systems (as each follow tighter relations at higher metallicities), but more likely it is indicative of other accretion events, whose debris may have been identified already in part (and discussed in e.g. Sec. 4.2.3). Because of the correlation between mass and metallicity, and between mass and star formation history, small mass objects will have lower average values of [Fe/H] (hence populate the [Fe/H] −1.3 regime), and also typically depict lower [α/Fe] due to their less efficient star formation. Both these facts seem to be consistent with the findings reported above. Nissen & Schuster (1997) were among the first to show that α-poor stars are interesting markers for accretion events, and as discussed in Ishigaki (2019), large surveys are helping to identify larger numbers of low-α stars which are then being followed-up with high-resolution (see e.g. Xing et al. 2019). Similarly also r-process enhanced stars are receiving more attention (see e.g. Sakari et al. 2018), particularly because the large scatter known to be present in the field halo population at very low metallicities has been suggested to be an indication of several accreted small galaxies (see e.g. Roederer et al. 2018). Although the formation channels of Carbon-enhanced metal-poor (CEMP) stars are not fully understood, they also provide interesting insights. For example, most CEMP-no stars (CEMP stars with no over-abundance of neutron capture elements) have [Fe/H] −2.5, while the CEMP-s stars (enhanced in s-process elements) typically have [Fe/H] > −2.5 (Aoki et al. 2007). Since s-process elements would be produced on longer timescales, such stars would have their origin (when not in a binary), in systems that have sustained star formation longer, i.e. more massive hosts, and thus be preferentially found in the inner halo (as the recent analysis of Lee et al. 2019, seems to support). On the other hand, CEMP-no would form more predominantly in low mass galaxies (hence their lower average metallicity), which when accreted would remain in the outer regions of the halo, because they would not be able to sink in via dynamical friction (see e.g. Starkenburg et al. 2017), which appears to be in line with the trends observed in the Milky Way halo (Carollo et al. 2014).
The general picture emerging, that of a few large building blocks dominating the inner stellar halo, the largest possibly being Gaia-Enceladus, is consistent with that expected from cosmological simulations for a galaxy like the Milky Way. Zoom-in cosmological hydrodynamical simulations such as AURIGA , FIRE (Hopkins et al. 2014), or those based on N-body methods combined with semi-analytic models such as Aquarius (Cooper et al. 2010), are very useful for understanding the general properties of stellar halos and how they might relate to e.g. merger history (Tissera et al. 2013, Pillepich et al. 2015, El-Badry et al. 2018. Some of these simulations are beginning to find look-alikes to the Milky Way also in terms of their merger history, as reported for example in Fattahi et al. (2019)/ Large fully cosmological hydrodynamical simulations such as the Illustris suite and its successor IllustrisTNG (Naiman et al. 2018) also contain Milky Waylike galaxies. The largest EAGLE cosmological simulation (Schaye et al. 2015) for example, contains 100 objects that are close analogs to the Milky Way in terms of stellar and dark mass, SFR and bulge-to-disk ratio, of which a handful with similar stellar halos in terms of their dynamical properties (see Bignone et al. 2019).
The evolution of the galaxy identified by Bignone et al. (2019) in the EAGLE suite to be a good Milky Way-alike is shown in Figure 16 around the time it merges with a Gaia-Enceladus analogue. The panels of the figure show the evolution of the star formation rate (top) and stellar mass (bottom) for the different "components" identified according to the circularity of their stars. At the time of the merger, there is a significant increase in the SFR in the whole system, with that of the fiducial thin disk increasing dramatically towards the end of the merger (fueled in part by gas from the accreted object, which also helps its further growth). Both panels show that all components are affected by the merger, suggesting the existence of populations in the thick disk and bulge that are co-eval with the timing of the merger.
The figure also shows that the Gaia-Enceladus-analogue barely completes two orbits before it is fully disrupted. Its mass ratio in this simulation is only 20%, and both the host as the infalling object have more than 50% of the baryons in cold gas. Such an event is thus quite different from those typically modeled in the context of dwarf galaxy accretion for two reasons. Firstly the debris will be much more complex because of the effect of dynamical friction: stars lost early on trailing or leading streams might have respectively, lower or higher eccentricities (and also lower metallicities). Their dynamical characteristics likely also differ from those lost in the next passage. Secondly, the gas will respond heavily, and  Evolution of a Milky Way-like galaxy identified in the largest EAGLE cosmological simulation. The galaxy's spheroid component has dynamical properties similar to those of the Milky Way (a significant population of stars on very eccentric orbits, i.e. the "Gaia-sausage") as a result of a merger with another system that was completed at z ∼ 1.2. Stars have been associated to "components" on the basis of their circularity (computed at each point in time), with low circularity (green) representing the spheroid, intermediate (orange) the thick disk and high circularity (blue) the thin disk. The thick disk in this simulation originates in part from a heated disk (star particles "transferred" from the thin disk and put on less circular orbits during the merger), but also from stars formed during the merger itself. Notice the increase in the SFR of all components during the merger. Credits: Figure  it is conceivable that some star formation might have taken place in the tidal arms or even the formation of star clusters to have been triggered. The degree of complexity of such an event evidences the need for tailored simulations including gas, star formation and chemical evolution to fully interpret and model the observations currently available.

Next steps: Simulations
Because of their limited resolution, many of the cosmological simulations just mentioned have been re-sampled to produce more particles (Lowing et al. 2015, Grand et al. 2018, Sanderson et al. 2018, as typically a star particle in a simulation nowadays should be seen as representing a stellar population of approximately 10 3 − 10 4 M . Nonetheless, there are still important limitations including the ability to trace the true phase-space distribution of stars originating in objects with a stellar mass lower than 10 6 M or thereabouts. Furthermore, not all these simulations follow chemical enrichment properly and this as discussed in previous sections, is necessary for the guidance in the interpretation of the various structures that are being (and will in the future, likely be) identified.
The study by Bignone et al. (2019) makes clear that cosmological simulations (even with low resolution) are useful also for exploring or making links between the formation paths of the different galactic components. Now that the identification of true analogues in cosmological simulations has become possible, it will be of great interest to carry out new zoom-in simulations of these objects. They will allow us to address a variety of questions, including for example gas physics, star formation and chemical enrichment processes, as well as to establish the true link between different events in the assembly history of the Milky Way. Furthermore, such simulations will be necessary to guide dark matter detection experiments which often assume that the dark matter particles follow Maxwellian velocity distributions. We now know that the stellar halo near the Sun is complex and has multiple kinematic components (see e.g. O'Hare et al. 2020, for a discussion of the impact on direct detection experiments). It has stars with kinematics corresponding to the tail of the thick disk, and from a component that is mildly retrograde which is associated to Gaia-Enceladus. But we do not know how the dark matter should be distributed given the particular Galactic history just unraveled (although see Necib et al. 2019a). It is therefore very important to carry out such simulations now that the amount of freedom has been significantly reduced and that the boundary conditions are better known so that we can provide concrete constraints on the initial conditions. Such simulations can be also used to understand some peculiarities about the Milky Way, such as the possibly fairly low mass for the super-massive black hole in the center of the Galaxy, or the distribution of satellites and the origin of their configuration.

Next steps: Statistical analyses
To trace the assembly history of the Milky Way as far back as possible, it will also be necessary to work in a more systematic fashion than until now. This is essential for establishing for example, the mass spectrum of the objects accreted and their internal characteristics. It will require the application and development of statistical methods, some of which are available in the literature (e.g. IoM, frequency space), Fourier analysis, machine learning, clustering algorithms such as (H)DBSCAN (as used by Borsato et al. 2020, Koppelman et al. 2019a, Necib et al. 2019b or based on neural networks (Yuan et al. 2019). An important aspect is the assessment of the statistical significance of a given feature or clump, and this can be done either through comparison to models or to suitably randomized samples. It would also be useful if there was more consistency in the different structures reported in the literature. Sometimes, a structure is reported as newly discovered, while it has been reported before (see O'Hare et al. 2020, for a recent example). This leads to confusion in the field (in the naming and in the reality of the features) and also does not help in building up a coherent picture of Galactic history. A possible improvement would be to assign membership probabilities to the different stars, and to publish these together with the different structures and stars' IDs.
Similar issues arise for the globular cluster population of the Milky Way. Although Massari et al. (2019) have tentatively associated many (at least 35%) of the globular clusters of the Milky Way to what we may identify as the main building blocks of the halo (at least near the Sun), namely Gaia-Enceladus, Sagittarius, the Helmi Streams and Sequioa, for many globular clusters the assignment is not unique, particularly in Integrals of Motion space (see for example Myeong et al. 2018b). More precise ages for the clusters could potentially aid as well because although the age-metallicity relations are well-defined, they are not fully unambiguous (see also Kruijssen et al. 2019).
For the interpretation of the various structures, it is clear that there is an urgent need for detailed chemical abundances of the stars with full phase-space information. Such samples will also aid in disentangling the different events, in assessing their origin (particularly if substructures overlap in e.g. IoM space, as will necessarily happen in the majority of cases), and characterizing their properties and history. A very nice example of what is currently possible along these lines was given in Das et al. (2019). Their analysis of APOGEE DR14 data is shown in Figure 17, where the different main structures discussed in this review are clearly separated in chemical abundance space.

Figure 17
Two-dimensional abundance distribution of APOGEE DR14 stars in different chemical abundance planes. The filled orange circles correspond to kinematically-selected disk stars from Bensby et al. (2014). Note the presence of Gaia-Enceladus (marked as "G-E" in all panels), and its clear distinction from the thick disk, not only in [Mg/Fe]

Next steps: Surveys
The need for spectroscopic surveys for large numbers of stars has long been recognized since the pioneering ideas of Freeman & Bland-Hawthorn (2002) and discussed in the context of e.g. ESO-ESA synergies by Turon et al. (2008). This interest has lead to a significant increase in the number and scope of the surveys. Surveys such as RAVE, SEGUE and Gaia-ESO 25 have been carried out over the past decade; some like APOGEE, LAMOST and GALAH have been running over the past years. In the next few years projects such as WEAVE, 4MOST and DESI using 4m class telescopes and MOONS 26 on the VLT will see the light. The Galactic surveys that will be carried out using these facilities have two complementary goals. The first is to obtain the missing radial velocity for stars for which Gaia has (at best) measured their 5D phase-space location. The second goal is to obtain high resolution spectra for brighter stars to do chemical labelling and characterization of the most metal-poor populations for example in the halo. The first goal is important for the dynamics of distant halo stars, for which tangential velocity constraints available from Gaia data releases will be available but less precise. The radial velocity measurements will allow mapping the mass distribution in our Galaxy at large radii . Also for nearby faint dwarf stars, obtaining the missing radial velocity component is useful because such stars could be used to trace the time-variations in the gravitational potential of the Milky Way for example, in frequency space as discussed in Sec. 3.3.
On the other hand, high resolution spectroscopic follow-up is of utmost importance as the discovery of Gaia-Enceladus has taught us. It is arguably the most powerful way to fully pin-down history and to identify debris with certainty, as well as to disentangle accretion events from one another. The high-resolution surveys planned by e.g. WEAVE and 4MOST will be carried out using relatively bright stars, i.e. down to G ∼ 16, because of the use of 4m class telescopes (see e.g. Christlieb et al. 2019). Although they will be invaluable, it is already clear that a wide-field spectroscopic survey on an 8-12m class telescope would be fantastic as it would really match the capabilities of Gaia, and of e.g. LSST 27 (Ivezić et al. 2019) in the coming years. In particular for the identification and characterization of debris from low mass systems, it will be necessary to target main sequence stars, as they are much more numerous than the few RGB stars present in ultra-faint-like galaxies. Such galaxies are particularly interesting because of questions related to the presence of thresholds for galaxy evolution, the impact of reionization and feedback processes, and because they may host some of the most metal-poor stars. These stars reveal the imprint of just a few supernovae and possibly of the initial mass function in the early Universe. Because of the low density of tidal debris and of the stellar halo more generally, follow-up must be carried out with wide-field, and an 8-10m telescope may well necessary to reach the required depth. The PSF 28 on Subaru (Tamura et al. 2016) is an instrument that could potentially help with the chemical labelling, although its highest resolution mode has R ∼ 5000 and so obtaining detailed chemical abundances for many elements will not be feasible. The MSE 29,30 is another interesting facility being considered, but there are no (other) concrete plans at the time of writing of this review, although Pasquini et al. (2018) discuss in some detail a concept developed at ESO whose main science driver is high-resolution follow-up of Gaia targets, in a case termed "the Milky way as a Model Galaxy Organism".

SUMMARY POINTS
1. The evidence discussed in this review seems to support the view that a milestone in Galactic history has been unveiled, namely the last big merger experienced by the Milky Way 10 Gyr ago. The merged galaxy, known as Gaia-Enceladus, appears to be responsible for a large fraction of the inner stellar halo. 2. A disk component was present at the time of merger with Gaia-Enceladus, which is possibly related to the metal-weak thick disk. This disk component was significantly perturbed as a consequence of the merger, and it is currently responsible for roughly 50% of the stars with halo-like kinematics near the Sun, i.e. it seems to constitute what has been known as the "in-situ" halo. There is also some evidence that the merger with Gaia-Enceladus triggered significant star formation in the thick disk, and it seems plausible that it was responsible for the formation of a large fraction of its stars. 3. These results appear to be consistent with the predictions of the ΛCDM model to first order. Objects in cosmological simulations with similar merger histories have been identified and are providing new insights in the evolution of the various Galactic components, and into how their histories may be interlinked. 4. Several other kinematic substructures have been identified amongst nearby stars likely associated to past accretion events. The Helmi streams (∼ 10 8 M ), (bonsai) Sequoia (∼ 10 7 M ), Thamnos (∼ 5 × 10 6 M ), together with Gaia-Enceladus (∼ 10 9 M ) and Sagittarius (∼ 5 × 10 8 M ), appear to be the largest building blocks of the halo of our Galaxy, although it is not always clear how and if they are related to each other. For some of these (Helmi streams, Gaia-Enceladus, and Sagittarius) it has been possible to derive (sketchy) star formation histories, and their chemical characterization has only just began via chemical labelling. The biggest current limitation is the availability of large samples of stars with detailed chemical abundance information, but when available this information turns out to be crucial. 5. When such samples become available it will be possible to do true Galactic archaeology, and establish the properties of the galaxies (star formation and chemical enrichment histories, mass, size) before they were accreted, i.e. really explore the high redshift universe from our own backyard. It should also be possible to carry out a similar exercise for the thick disk, or the proto-disk, which we know now was present 10 Gyr ago (at z ∼ 1.8), and seems to have been traced back to metallicities [Fe/H] −4 (Sestito et al. 2019). Establishing the properties of this disk (structure, star formation) will allow us to make a direct link to high-redshift disks being currently observed in-situ. 6. Most of the above mentioned building blocks were identified using full phase-space information of stars in the Solar vicinity. At larger distances (beyond 20 kpc from the Galactic center, 10 kpc from the Sun), little is known about kinematics and chemistry of the stars but several large substructures have been known for a while thanks to wide-field photometric surveys. These include Hercules-Aquila and the Virgo overdensity, as well as several others. The relation between these and the so-far identified building blocks if any, has not been really pinned-down although we present here a way forward using orbital integrations. With Gaia DR3 and subsequent releases (DR4 and beyond, since the mission has been extended over the nominal lifetime for at least 2 years), and in combination with spectroscopic surveys, it should be possible to address this and many more questions. 7. Also in the direction towards the Galactic center, little is known, yet this is where the halo reaches its highest density. Merger debris from another massive building block might well be present (as suggested by the analysis of the age-metallicity relations and dynamics of globular clusters by Kruijssen et al. 2019, Massari et al. 2019, and which could be a "Kraken"-like object). It is not clear at this point if Gaia can answer this, but an astrometric mission in the near infrared (see e.g. Gouda 2015, Hobbs et al. 2019) could perhaps address this open question. 8. There is more to be gained from the analysis of substructures in the vicinity of the Sun, and some may well be due to internal mechanisms (such as the Galactic bar or non-integrability of a generic Galactic potential), or even reveal the response of the Galaxy to an accretion event. Such analysis as well as addressing the above mentioned items, will require more detailed modeling, also from a dynamical perspective, and preferably through zoom-in cosmological simulations which now could model systems which at least in terms of their merger history, would be much more representative of the Milky Way.
Enormous progress has been made in recent years in our understanding of the evolution of the Milky Way from the perspective of the stellar halo and thick disk disk. This has been possible particularly thanks to Gaia DR2 and the many spectroscopic surveys currently available, and especially those surveys with detailed chemical abundance information have been proven to be crucial. Yet a plethora of questions remain open for the next decade. Fortunately we already know we will have many of the necessary tools to answer them. There are very exciting times ahead of us in the field of Galactic archaeology.

FUTURE ISSUES
1. Although challenging, it would be extremely interesting to identify accretion events that have taken place even before the merger with Gaia-Enceladus 10 Gyr ago. 2. The characterization (in terms of dynamical, chemical and star formation history) of the disk present at the time of the merger with Gaia-Enceladus, should be pursued. An important link is there to be made between the disks observed in-situ in high-z observational studies, and that revealed in the Milky Way. 3. Zoom-in cosmological simulations of systems with a merger history and dynamics similar to that suggested by recent data would be particularly useful for understanding Galactic history and the link between its various components and detailed properties. Such simulations would also allow to make robust predictions for direct detection dark matter experiments. 4. There is an ever increasing need for large high-resolution spectroscopic surveys of relatively faint stars (G 16) to supplement the dynamical information that is becoming available thanks to the Gaia mission. Determination of precise ages for such a sample of stars would be highly valuable and useful to date various events in Galactic history, particularly at early times.

DISCLOSURE STATEMENT
The author is not aware of any affiliations, memberships, funding, or financial holdings that might be perceived as affecting the objectivity of this review.