© ESO, 2017

1. Introduction

Current spectroscopic surveys, such as APOGEE, APOGEE-2, GES, GALAH, and RAVE (Allende Prieto et al. 2008; Majewski et al. 2010; Eisenstein et al. 2011; Sobeck et al. 2014; Gilmore et al. 2012; Anguiano et al. 2014; De Silva et al. 2015; Steinmetz et al. 2006; Zwitter et al. 2008; Siebert et al. 2011; Kordopatis et al. 2013), are extending our knowledge of the chemical composition and radial velocities of several hundred thousand stars in the Galaxy up to several kpc from the Sun. In the coming years, starting from the first release in mid-September 2016, the ESA astrometric mission Gaia (Perryman et al. 2001) will deliver parallaxes and proper motions for approximately 1 billion stars in the Galaxy and radial velocities for approximately one tenth of them. This unprecedented amount of data, coupled with Gaia follow-up spectroscopic surveys currently under definition (WEAVE, 4MOST, and MOONS (Dalton et al. 2012; de Jong et al. 2012; Cirasuolo et al. 2012) will potentially allow us to answer many simple, but still unanswered, questions concerning Galactic studies and, more specifically, topics such as the characteristics of the different Milky Way stellar populations, the mechanisms by which they have been shaped over time, their evolutionary links, and ultimately the extent to which they are the result of in-situ star formation or if they are simply the deposit of stars and structures accreted over time. The Milky Way will constitute a natural testbed for any cosmological scenario, in particular for Λ cold dark matter (CDM) cosmology, according to which galaxies grow hierarchically, from the formation of small subunits that subsequently merge to evolve into the galaxies that we see today (White & Rees 1978).

ΛCDM models predict indeed that a galaxy such as the Milky Way should contain hundreds of stellar streams in the solar vicinity (Helmi et al. 1999, 2003; Gómez et al. 2013), relics of the merging of other galactic systems, with masses comparable to or significantly smaller than our own Galaxy at the time of their accretion (Bullock & Johnston 2005; Stewart et al. 2008; De Lucia & Helmi 2008; Cooper et al. 2010; Font et al. 2011; Brook et al. 2012; Pillepich et al. 2015; Rodriguez-Gomez et al. 2016; Deason et al. 2016). While we have evidence of recent and ongoing accretion events experienced by our Galaxy (Ibata et al. 1994, 2001; Belokurov et al. 2007) and some streams – vestige of more ancient mergers – have been possibly detected also in the solar vicinity (Helmi et al. 1999, 2006; Klement et al. 2009; Smith et al. 2009; Nissen & Schuster 2010), we are still far from the numbers predicted by cosmological simulations. This may be an indication that we have overestimated the importance of accretions in building a Galaxy like our own over time, or this may be rather a consequence of the difficulty in finding the remnants of accreted satellites, in particular of those related to the most ancient accretion events, once they are spatially mixed with in-situ stars in the Milky Way.

Numerous studies in the last 15 yr have suggested that the imprints of past merger events, even if dispersed in configuration space, can still be identified in kinematics-related spaces, such as the energy-angular momentum space (hereafter referred to as E − L_z, L_z being the z component of the angular momentum for an axisymmetric potential where z is the symmetry axis), the L_⊥ − L_z space, where L_⊥ is the angular momentum component in the x − y plane, or spaces of orbital parameters, such as the apocenter-pericenter-angular momentum space (APL) and its projections (Helmi et al. 1999, 2006; Helmi & de Zeeuw 2000; Knebe et al. 2005; Brown et al. 2005; Font et al. 2006; Choi et al. 2007; Morrison et al. 2009; Gómez et al. 2010; Re Fiorentin et al. 2015) – see also Klement (2010), Smith (2016) for two recent and comprehensive reviews.

For example, distinct accreted satellites should appear as coherent structures in the E − L_z space (Helmi & de Zeeuw 2000) and the shape of these structures should not significantly change during the accretion event, even in the case of a time-dependent potential (Knebe et al. 2005; Font et al. 2006; Gómez et al. 2010). From the number of clumps found in the E − L_z space of several million Gaia stars, it should potentially be possible to set a lower limit to the number of past accretion events (Helmi & de Zeeuw 2000; Gómez et al. 2010). But it is still necessary to check whether this search is really feasible and meaningful, and what an overdensity in the E − L_z space really means. Can we realistically make use of the number of substructures found to set a lower limit to the number of merger events that the Galaxy experienced in its past ?

Other spaces where the signature of accretion events may be found are the apocenter-pericenter space: here accreted stars should tend to cluster around a common value of eccentricity (Helmi et al. 2006), even if they can span different apocenter and pericenter distances. Here we are interested in understanding if this coherence in the AP space is really maintained, and again, if it is possible to separate accreted stars from the underlying, in-situ population, which – as we show – can be also clustered in this space.

Finally, about the L_⊥ − L_z space, we wonder what are the regions in this space where the probability of finding accreted stars is the highest.

In this paper, we address these questions, by means of high-resolution N-body simulations of a Milky Way-type galaxy which accretes one or several (up to four) satellites. The novelty of this work with respect to the majority of those already published, is our analysis of fully consistent models, where both the satellite(s) and the Milky Way galaxy are “live” systems, which can react to the interaction, experience kinematical heating, tidal effects and dynamical friction (the latter, a process often neglected in previous studies, but see Knebe et al. 2005; Meza et al. 2005; Font et al. 2006; Re Fiorentin et al. 2015), and exchange energy and angular momentum. We show how, differently from previous findings, the energy-angular momentum space – and similarly, the L_⊥ − L_z and the apocenter-pericenter-angular momentum spaces – become hardly decipherable spaces when some of the most limiting assumptions, which have affected previous works, are removed (see Sect. 2.1). We show indeed that:

• because the energy and angular momentum of a satellite are notconserved quantities during an interaction, each satellite givesorigin to several independent overdensities;

• multiple satellites overlap in the E − L_z space;

• in-situ stars affected by the interaction(s) tend to progressively populate a region of lower angular momentum and higher energy that the one initially (i.e., before the interaction) populated;

• most of the accreted stars overlap with in-situ stars.

We conclude from these findings that the search for the building blocks of our Galaxy in kinematic spaces is mostly inefficient. Similar conclusions are found for the APL and L_⊥ − L_z spaces.

In agreement with previous works (see also Gould 2003), we find that all these spaces are rich in substructures, but that the origin of these substructures cannot be determined with kinematics alone. As an example, we show here that the in-situ stellar halo, formed as a result of the interaction, is neither smooth nor non-rotating. This is different from what is usually assumed for in-situ halo stars in the Galaxy (Helmi et al. 1999; Kepley et al. 2007; Smith et al. 2009). Moreover, in-situ halo stars can also be clustered in kinematic spaces, and extreme caution must be employed before interpreting over-densities in those spaces as evidence of relics of past accretion events (see, for example, Morrison et al. 2009; Helmi et al. 2006; Klement et al. 2009; Re Fiorentin et al. 2015). Detailed chemical abundances and/or ages will definitely be necessary to identify stellar streams that should show up in those spaces as patterns distinct from those described by the in-situ stellar populations.

This paper is organized as follows: in Sect. 2 we present the N-body simulations analyzed in this paper (Sect. 2.2) after discussing the main numerical works published on the subject so far and the approaches previously adopted (Sect. 2.1). In Sect. 3, we present the main results of our work, starting by analyzing the E − L_z space (Sect. 3.1) and then move on to show the results for the L_⊥ − L_z (Sect. 3.2) and the APL (Sect. 3.3) spaces. A discussion of our findings is provided in Sect. 4, and finally, in Sect. 5, the main conclusions of our work are presented.

2. Numerical methods

2.1. Previous numerical modeling

Before describing the N-body simulations analyzed in this paper, it is worth summarizing the main characteristics and the main dynamical processes modeled so far in the literature to study the kinematic signatures of accreted stars in the Galaxy. As we see indeed highlight here, understanding these characteristics is essential to understanding the reasons why the present work does not reach the same conclusions reached by previous studies.

Schematically, the suggestion to detect stellar streams in the Milky Way by looking for (sub)-structures in kinematic spaces is mainly based on two kinds of models; test particle models and self-consistent, N-body models.

Test particle models have been developed since the pioneering works of Helmi et al. (1999), Helmi & de Zeeuw (2000). In these studies (see also Kepley et al. 2007) the Milky Way is represented by a fixed, rigid potential and the satellites as a collection of particles, usually gravitationally interacting. In most recent works, the Milky Way potential has been allowed to vary with time (see, for example Gómez & Helmi 2010; Gómez et al. 2010). All these models have confirmed Helmi’s early suggestions: while the debris of early accreted satellites is very difficult to recover spatially, strong correlations and structures should be visible in velocity spaces (E − L_z, L_⊥ − L_z, and APL). We note that because these models assume an analytic Galactic potential and no dynamical friction is taken into account, the energies and angular momenta of the accreted satellites will be overall conserved by definition (for the angular momenta, L, this is true for a spherical potential, while in an axisymmetric one, it will be only the z −component to be strictly conserved.). Thus the finding that accreted satellites stand out as lumps in integrals of motion spaces and in particular that each lump corresponds to an accreted satellite (see for example Helmi & de Zeeuw 2000, for the E − L_z space), is, generally, a direct consequence of the assumptions made in these models, rather than an intrinsic feature of the accretion event. Unless satellites disperse before suffering substantial dynamical friction – but note that in this case they would hardly attain the inner Galactic halo (distances from the Galaxy center ~15 − 20 kpc, Carollo et al. 2007), where kinematic detections of streams are currently mainly focused – satellites lose their coherence in integrals of motion spaces during their accretion into the Galaxy (see Sect. 3).

Caution should also be taken due to the fact that these models generally do not include any in-situ (i.e., formed in the Galaxy) population, that acts, as we see here, as a background where the signal of accretion events is mostly washed out. In all the above cited works, indeed, the stellar halo is exclusively made of accreted stars. This assumption has the natural consequence of maximizing their signal, even when observational errors are taken into account, as in Helmi & de Zeeuw (2000) and Gómez et al. (2010). An exception to this general approach is represented by the work by Brown et al. (2005), who also modeled an in-situ population, concluding that a search for clumpy structures in the E − L_z space is indeed very challenging for astrometric surveys like Gaia, one of the reasons being the presence of the background Galactic population. We note, however, that Brown et al. (2005) modeled the in-situ halo population as a smooth distribution in phase-space, added a posteriori to the model: as a consequence, this population cannot respond self-consistently to the interaction, or experience either kinematical heating or any clustering in kinematic spaces. The natural consequence is that the lumpiness of the accreted population is still overestimated, because it is modeled against an assumed smooth background. We show in the following pages that the situation is even more complicated than that discussed by Brown et al. (2005) and the search even more inefficient, when, also, the in-situ population is modeled self-consistently.

Self-consistent, “live” N-body models remove the hypothesis of a rigid Galactic potential, since both the Galaxy and the satellite(s) are modeled as a collection of particles (dark matter and/or stars) that respond self-consistently to the merger. Dynamical friction is thus taken naturally into account and so should also be the response of the in-situ Galactic populations to the accretion events. To our knowledge, the importance of using this kind of approach for detecting streams kinematically was first pointed out by Knebe et al. (2005). They used cosmological, dark matter-only simulations and showed that the lumpy appearance of accreted satellites is significantly smeared out when a “live” model is adopted. Unfortunately, Knebe et al. (2005) do not discuss nor show the response of the in-situ component to the interaction and, as a consequence, the efficiency of the “integrals-of-motion” approach is still overestimated in their work. Many other simulations have investigated the E − L_z space, or equivalent, by making use of self-consistent, N-body simulations (see, for example, Meza et al. 2005; Helmi et al. 2006; Font et al. 2006; Gómez et al. 2013; Re Fiorentin et al. 2015). However, either they have not discussed the response of the in-situ population to the accretion event(s) (Meza et al. 2005; Helmi et al. 2006; Font et al. 2006; Gómez et al. 2013), or they have not taken into account the in-situ halo population that naturally forms in merger events (as shown, for example, by Zolotov et al. 2009; Purcell et al. 2010; Font et al. 2011; Qu et al. 2011; McCarthy et al. 2012; Cooper et al. 2015), replacing it with a smooth halo component added a posteriori (Re Fiorentin et al. 2015). To our knowledge, the only works that have started to investigate the question of the overlap between in-situ and accreted stars are those by Ruchti et al. (2014, 2015), however they investigated mostly low-mass mergers (mass ratio ~1:100).

The limitations that, in our opinion, nearly all previous works have suffered constitute the main motivation of our work. We consider it to be the first attempt to model both the accreted and in-situ stellar populations during one or multiple accretion events in a more realistic way. This is an unavoidable step in making predictions about the redistribution of in-situ and accreted stars in the Galaxy in view of Gaia, but also in interpretation of current kinematic data where substructures are currently detected and an extra-galactic origin is often favored.

2.2. Our simulations

To this end, published here are the analyses of three dissipationless, high-resolution simulations of a Milky Way-type galaxy accreting one, two or four satellites, respectively. By also including a 4 × (1:10) merger in the present paper we do not mean to suggest that the Milky Way should contain such a high fraction of ex-situ stellar mass but we consider that the analysis of this simulation brings an important element to the discussion and to the comprehension of the overlap of in-situ/ accreted stars. One may naively think that accreting more mass can increase the signal of accreted over in-situ stars in kinematic spaces. By adding an extreme case, such as that of the 4 × (1:10) merger, we show in the following sections that this is not the case: it is true that the remnant galaxy ends up with a larger fraction of accreted compared to in-situ material, but at the same time (1) this material is redistributed over a larger portion of the phase-space and (2) the amount of in-situ stars heated by the interaction and found in the halo region also increases. The final result is that even a very extreme situation such as that of a 4 × (1:10) merger would not be more favorable to the kinematic detection of streams than a simple 1:10 or 2 × (1:10) merger. We consider this an important element of the following discussion.

Each satellite has a mass which is one tenth of the mass of the Milky Way-like galaxy. The total number of particles used in these simulations varies between N_tot = 27500110, for the case of a single accretion, and N_tot = 35000140, when four accretions are modeled. The massive galaxy contains a thin, an intermediate, and a thick stellar disc (mimicking the Galactic thin disc, the young thick disc and the old thick disc, respectively) (see Haywood et al. 2013; Di Matteo 2016), a population of a hundred disc globular clusters represented as point masses, and is embedded in a dark matter halo. The total number of stellar (thin and thick disc) particles in the main galaxy is 20000000, the number of globular cluster particles is 100 and the number of dark matter particles is 5000000. The discs are modeled with Miyamoto-Nagai density distributions, of the form

(1)

with M_∗,a_∗ and h_∗ masses, characteristic lengths, and heights which vary for the thin, intermediate, and thick disc populations (see Table 1 for all the values); the system of disc globular clusters has scale length and scale height equal to those of the thick disc and the dark matter halo is modeled as a Plummer sphere, whose density is: (2)with M_halo and a_halo characteristic mass and radius, respectively (see Table 1 for all the values). The satellite galaxies are rescaled versions of the main galaxy, with masses and total number of particles ten times smaller, and sizes reduced by a factor (again, see Table 1). It has been shown that most of the initial dark halo mass is rapidly lost by the satellite galaxy as it is placed in the host potential; Villalobos & Helmi (2008), for example, estimate that after a few 10⁸ Myr, that is, before the first pericenter passage, only ~30% of the initial dark halo mass is still bound to the system. We have assumed an initial halo mass for the satellite of M_halo = 1.6 × 10¹⁰M_⊙ and a visible mass of M_⋆ = 5 × 10⁹M_⊙, which implies a M_⋆/M_halo ratio = 0.3, compatible with the M_⋆/M_halo ratio ~0.2 found in the simulations by Villalobos & Helmi (2008) before the first close encounter between the satellite and their modeled Milky Way-type galaxy. We also highlight that our choice of a core dark matter halo for both the Milky Way-type galaxy and the satellite(s) was prompted by a number of observations that seem to be more consistent with a dark halo profile with a nearly flat density core (Flores & Primack 1994; de Blok et al. 2001b,a; Marchesini et al. 2002; Gentile et al. 2005; Kuzio de Naray et al. 2006, 2008).

The Milky Way disc initially coincides with the x − y plane of the reference frame with the spin of the Milky Way-type galaxy being oriented as the z-axis and positive. In a reference frame with the origin at the center of the Milky Way-type galaxy, the x − y plane coinciding with its disc and the z-axis being oriented as the spin of the Milky Way-type galaxy, satellite positions are (x_sat,y_sat,z_sat) and their distances from the center of the Milky Way-type galaxy are equal to D_sat. Their initial 3D velocities, relative to the Milky Way-type galaxy are (v_x,sat,v_y,sat,v_z,sat), V_sat being their absolute value. Their orbital planes are inclined by (θ_orb,φ_orb) with respect to the Milky Way-type galaxy, θ_orb being the angle between the intersection of the satellite orbital plane with the xy-plane and the x-axis and φ_orb being the angle between the projection of the satellite orbital plane onto the xz-plane and the x-axis. Each satellite has an internal angular momentum whose orientation in space is described by the angles (θ_sat and φ_sat) with θ_sat being the angle between the projection of the satellite spin onto the xy-plane and the x-axis and φ_sat being the angle between the projection of the satellite spin onto the xz-plane and the x-axis. A schematic representation of all these angles is given in Fig. 1. Their values are reported in Table 2 together with positions, velocities, and initial orbital angular momenta. We note, in particular, that:

• all satellites are initially on direct orbits(L_z,sat, the z −component of their orbital angular momentum, being positive), except for satellite 4 of the 4 × (1:10) merger, which is initially placed on a retrograde orbit, L_z,sat being negative (see Table 2).

• Initial orbital velocities of the satellite galaxy correspond to those of a parabolic orbit for a 1:10 mass ratio with a MW-type galaxy mass of 100 (in our mass units). This value for the total mass is slightly higher (~10%) than the initial total mass in the Milky Way-type galaxy but since the Milky Way-type system in our models accretes not only one, but 2 or 4 satellites, we have taken a slightly higher initial mass to also take into account (at least partially) this mass growth. The choice of exploring parabolic orbits, in particular, is in agreement with cosmological predictions (Khochfar & Burkert 2006).

We comment further on the choice of our initial conditions in Sect. 4.5.

Initial conditions were generated adopting the iterative method described in Rodionov et al. (2009). All simulations were run by making use of the TreeSPH code described in Semelin & Combes (2002), which has already been extensively used by our group over the last ten years to simulate the secular evolution of galaxies and accretions and mergers as well. Gravitational forces are calculated using a tolerance parameter θ = 0.7 and include terms up to the quadrupole order in the multiple expansion. A Plummer potential is used to soften gravitational forces, with constant softening lengths for different species of particles. In all the simulations described here, we adopt ϵ = 50 pc. The equations of motion are integrated using a leapfrog algorithm with a fixed time step of Δt = 2.5 × 10⁵ yr.

In this work, we make use of the following set of units: distances are given in kpc, masses in units of 2.3 × 10⁹M_⊙, velocities in units of 100 km s^-1 and G = 1. Energies are thus given in units of 10⁴ km²/ s² and time is in units of 10⁷ yr. With this choice of units, the stellar mass of the Milky Way-type galaxy at the beginning of the simulation is 5.1 × 10¹⁰M_⊙. Since all the N-body models presented in this work are dissipationless, units can easily be re-scaled and, for example, masses can be reduced to mimic a Milky Way-type galaxy at higher redshifts. Such a rescaling would still be within estimated uncertainties, even for dark matter halos with masses below 10¹²M_⊙, that is in the steep part of the M_⋆/M_halo relation, at least between z = 0 and z = 1. Indeed, for a galaxy with a halo mass M_halo ~ 2 × 10¹¹ at z = 0, and a corresponding M_⋆/M_halo ~ 0.01 at that redshift, if one decreases its dark matter mass by a factor of 2 between z = 0 and z = 1 and at the same time maintains the M_⋆/M_halo ratio constant and equal to 0.01, this ratio would be between 2 and 3 sigma from the mean relation at z = 1 (see Figs. 7 and 8 from Behroozi et al. 2013).

3. Results

In the following analysis, unless explicitly stated, all quantities are evaluated in a reference frame whose origin is at the centre of the Milky Way-type galaxy. This center is evaluated, at each snapshot of the simulations, as the density center, following the method described in Casertano & Hut (1985). The x − y plane coincides with its disc and the z-axis is perpendicular to it. All the results presented in this section concern the stellar distribution in integrals-of-motion and kinematic spaces, in the case of a Milky Way-type galaxy accreting one or several satellites. We refer the reader to Appendix A for the distributions obtained when the Milky Way is evolved in isolation without experiencing any accretion.

3.1. On the E – L_z space

We begin our investigation by looking at the E − L_z space, which has been proposed as a natural space where to look for the signatures of past accretion events (Helmi & de Zeeuw 2000). The results are structured in the following way: first we discuss how one or several satellites redistribute their stars in this space during their accretion into the galactic potential, then we show the predictions of our models concerning the overlap between accreted and in-situ stars, and finally we restrict our analysis to the solar vicinity and discuss the difficulty and overall efficiency of searching for stellar streams in the E − L_z space.

3.1.1. Coherence of accreted structures

The redistribution in the E − L_z space of satellite stars during their accretion into our Milky Way-type galaxy is shown in Fig. 2 at different times during a 1:10 merger. The satellite represents, at the beginning of the interaction, a clump in the E − L_z space, characterized by high energy and a large spread in the z −component of the angular momentum. As it approaches the first pericenter passage, at t = 0.5 Gyr, the distribution of satellite stars in the E − L_z space changes: stars are redistributed over a much larger interval of energies, from slightly positive (for stars which become unbound to the system), to significantly negative (for those stars which approach the innermost galactic regions). At the same time, the satellite becomes closer overall to the galactic center and dynamical friction becomes stronger: as a consequence, the angular momentum tends to diminish, hence the increased compact distribution in L_z found at this epoch. This strong elongation in the redistribution of stars in the E − L_z space is reminiscent of the strong spatial elongation of satellite galaxies (and even less massive systems, such as globular clusters) when they approach the pericenter of their orbit: tidal tails are particularly elongated and relatively thin at this epoch and can span up to tens of kpc around the satellite bound mass. As the satellite reaches its next apocenter passage, its spatial distribution becomes more compact, as does its energy distribution (see Fig. 2 a time t = 1 Gyr), while the z −component of its angular momentum becomes broader again¹. At each orbit, satellite stars go through this phase of compression and expansion in the E − L_z space. Globally, as an effect of dynamical friction, the energy and the angular momentum decrease and the satellite penetrates deeper and deeper the potential well of the main galaxy. During the interaction, a portion of the stars becomes unbound, leaves the satellite and goes to populate the tidal tails which develop around it. Because stars lost at different passages are characterized by different values of energy and angular momenta (at first approximation, their energy and angular momenta at the moment they escape the satellite are those of the satellite center), escaped stars pack into different substructures, depending on the time they escaped the satellite. Stars lost at early epochs populate the upper part of the E − L_z diagram, while stars lost at more advanced epochs of the interaction are preferentially grouped into substructures of lower (i.e., more bound) energies. Thus, in general, if dynamical friction has time to act on the satellite before it becomes a gravitational unbound set of stars, satellite stars loose their coherence in the E − L_z space (as also found by Meza et al. 2005): a satellite gives rise to several streams or horizontal ripples, which are not homogenous all along their length and can become clumpy in the middle. In the following, we use the terms “ripples” or “streams” to identify the large-scale over-densities in the E − L_z plane and refer to “clumps” to identify either the high-density regions of those streams, or small-scale over-densities. The number and density of these ripples depend on the number of passages the satellite experienced around the main galaxy and on the mass-loss it experienced at each passage. Note that the loss of coherence of satellite stars in E − L_z space is not dependent on the particular choice of orbital parameters: in all simulated cases, from the case of a single 1:10 accretion, to the case of 2 × (1:10) and 4 × (1:10) mergers, each satellite contributes to several ripples and clumps in the E − L_z space, redistributing its stars to a large extent, both in energy and angular momentum (see Fig. 3). Note also that when multiple satellites are accreted, even if their initial energies and/or angular momenta are different, once the merger is completed their stars tend to redistribute over a similar portion of the E − L_z space, that is, the overlap between stars initially belonging to different satellites is not negligible. This indicates that, even in the ideal case of the absence of an in-situ population, a given region of the E − L_z space and in particular any given stream in this space, can be the result of the overlap of different accreted structures. Finally, it is worth emphasizing that accreted stars do not only redistribute in clumps: part of them is more diffusely distributed and cannot be associated with any clear over-density (see Fig. 3).

3.1.2. Overlap with in-situ stars

In the previous part of this section we saw how one or several satellites accreted onto a Milky Way-type galaxy redistribute in the E − L_z space. We saw that ripples and clumps in this space cannot be associated to a single accretion event: a satellite, during its accretion into the Galaxy, redistributes in several streams, and different satellites, with different initial orbital parameters, show a not negligible overlap in this space. However, there is another element that complicates the research of substructures in the E − L_z space: the presence of in-situ stars that respond to the interaction and redistribute in the energy-angular momentum plane, occupying a region similar to that of accreted stars, as we discuss further below.

It has been known for two decades that in-situ disc stars, kinematically heated during accretion events, redistribute in a thicker disc and inner halo (Quinn et al. 1993; Walker et al. 1996; Villalobos & Helmi 2008, 2009; Zolotov et al. 2009; Purcell et al. 2010; Font et al. 2011; Di Matteo et al. 2011; Qu et al. 2011; McCarthy et al. 2012; Cooper et al. 2015). As a result of this kinematic heating, in-situ stars tend to gain some energy and lose part of their initial angular momentum, since their orbits tend to become more radially and vertically elongated. This change in the spatial and kinematic properties of in-situ disc stars during an interaction also has some consequences, naturally, for their redistribution in the E − L_z space, as shown in the central column of Fig. 4. Stars that, in the absence of an interaction, would remain confined to a relatively thin and elongated region of the E − L_z space corresponding to that occupied by our initial thin/thick disc (see Fig. A.1), tend to redistribute toward lower levels of angular momentum and higher energies (see Fig. 4). The higher the number of accreted satellites (and thus the larger the accreted mass), the broader the distribution of in-situ stars in E − L_z space (Fig. 4). Moreover, the distribution of in-situ stars is clumpy in this space. In the case of a single accretion, as well as in the case of the accretion of four satellites, clumps appear not only in the region occupied by the stellar disc, as already pointed out by Gómez et al. (2012), but also in the less bound region of the diagram – the region naturally occupied by halo stars (an example of the spatial distribution of stars belonging to some of these clumps, for the case of the accretion of two satellites, is given and discussed in Appendix B). When all stars are taken into account, without any differentiation on their origin (in-situ or accreted), one sees clearly that this space is barely decipherable: part of the clumps found (see e.g., those in the bottom right-hand panel of Fig. 4, at 0 <L_z ≤ 20 and − 5 ≤ E< 0) do not have an extragalactic origin, but are made of in-situ stars. This finding is valid globally, but also locally, that is, when the search is restricted to solar vicinity volumes, as we show in the following.

Before moving to local searches for accreted streams, however, we want to direct the reader’s attention to another point that seems to us worth emphasizing: as a consequence of the angular momentum redistribution taking place during mergers, some of the in-situ disc stars, that initially rotate in a direct sense, with a positive L_z, tend to redistribute onto retrograde orbits. The larger the accreted mass, the higher the number of in-situ stars that end up with negative L_z and the greater the maximum value of the negative angular momentum attained. In the case of 4 × (1:10) mergers, for example, some high-energy, in-situ stars are on retrograde orbits with angular momenta lower than L_z = − 20 (in our units, see Sect. 2.2). As a consequence, we caution that the presence of stars on retrograde orbits in the Galaxy, even at high energies, is not necessarily evidence of an accreted, extragalactic origin.

3.1.3. Looking for streams in the solar vicinity

After discussing the global behavior of accreted and in-situ stars in the E − L_z space, we now investigate the feasibility of detecting stellar streams in this space, in local volumes. Defining “solar vicinity” volumes in an N-body simulation is not trivial. One can choose to place the local volume either at a distance from the galaxy center comparable to the Sun–Galactic center distance (i.e., 8–8.5 kpc), or to relate its position to the length of the stellar bar, or to the bar’s resonances. Moreover, because accreted stars do not necessarily redistribute homogeneously in configuration space, the fraction of accreted stars in a particular solar volume can significantly vary depending on the particular choice of that volume. To try to make our analysis as generic as possible and independent from the particular choice of our solar vicinity volume, we have chosen to identify several regions in our simulated galaxies, as schematically represented in Fig. 5: for each simulation, we have indeed defined 16 spherical volumes, all with radii of 3 kpc, centered at a distance of 8 kpc or 12 kpc from the galaxy center and distributed homogeneously in azimuth. For all stars in each of these volumes at the final time of the simulation, we have analyzed their distribution in the E − L_z space. For two of these volumes, one placed at 8 kpc and one at 12 kpc from the galaxy center, the corresponding distribution is shown in Fig. 6, for the case of a single merger (i.e., 1 × (1:10)). Despite some differences in the values attained by the stellar distribution in the different regions one notes, for example, that in the volume at 12 kpc, for any given value of energy, stars attain larger values of L_z, as expected, some points are common to all examined cases:

• 1.

  in a 3 kpc–wide region around the Sun, theE − L_z space is rich in substructures;

• 2.

  these substructures are present in the in-situ population and in the accreted population as well;

• 3.

  a single satellite gives rise to a multitude of substructures, of different extents and sizes;

• 4.

  the overlap between the in-situ population and the accreted one is substantial, to the point that no clear or evident distinction can be made between the two on the basis of the analysis of the E − L_z space alone;

• 5.

  in none of these volumes, can the presence of a substructure be easily and unambigously interpreted as evidence of an extragalactic origin for the stars that compose it: several stellar clumps, also some of those at moderate or negative L_z and relatively high energies, have indeed an in-situ origin, being made of stars initially in the disc, then heated by the interaction.

The analysis presented in Fig. 6 has been made considering all stars (i.e., disc(s), as well as halo stars) in a given stellar volume. One may thus think to alleviate the problem of the contamination of the in-situ component by kinematically selecting halo stars only. This is indeed a common strategy, supported by the idea that it is among the oldest populations in our Galaxy – currently found in the stellar halo – that the remnants of past accretion events can be found. We have thus repeated the analysis presented in Fig. 6, this time by only applying it to kinematically defined halo stars. To this end, for each solar vicinity volume we have selected only stars whose velocities (U_LSR,V_LSR,W_LSR) in the Local Standard of Rest (i.e., LSR) reference frame satisfy the condition km s^-1. Here U_LSR,V_LSR and W_LSR denote, respectively, the radial, tangential and vertical velocity in a reference frame whose origin is in the chosen stellar particle position, R_⋆ and which rotates around the galaxy center at a constant circular velocity, V_Sun, equal to V_C(R = R_Sun), i.e. equal to the value of the rotation curve of the modeled galaxy at R = R_Sun. The threshold of 180 km s^-1 is commonly used to select halo stars in observational samples at the solar vicinity (see, for example Nissen & Schuster 2010). Moreover, because in the selected “solar vicinity” volumes, our modeled galaxies attain values of the rotation curve similar to those estimated in our Galaxy at the solar radius (i.e. between 200 and 250 km s^-1), we feel comfortable in applying the same selection criterion also to our modeled data. This selection is shown in the top left-hand panel of Fig. C.1, for the same volume considered in the upper row of Fig. 6. As it can be appreciated from this figure that even when the analysis is restricted to kinematically selected halo stars only, the two main problems affecting the analysis presented in Fig. 6 are not alleviated:

• 1.

  the contamination of in-situ stars is still severe;

• 2.

  the clumpy redistribution of stars in the E − L_z space is not a characteristic of the accreted population only, since in-situ halo stars also show an inhomogeneous distribution.

In fact, accreted stars dominate only a marginal region of the E − L_z space of halo stars (see bottom-middle panel of Fig. C.1). Most of this space, which is all but smoothly distributed, is indeed dominated by the in-situ halo population (see bottom-right panel of Fig. C.1). To better quantify the importance of accreted stars, we have evaluated the area of the region, in the E − L_z space, occupied by accreted stars. To this end, we have selected only regions where the percentage of accreted stars is greater or equal to 90% of the total (this corresponds, for example, to the red-brown regions in Fig. C.1, bottom-middle panel) and we have divided this area by the area of the E − L_z space occupied by halo stars. We refer to this fraction as f₉₀. For the volume shown in Fig. C.1, f₉₀ = 0.16. Thus, in this volume, accreted stars dominate only a marginal region of the overall distribution. We can relax the threshold and look for the area where the percentage of accreted stars is greater or equal to 60% of the total distribution and, as before, divide this area by that occupied by the whole distribution. We call this fraction f₆₀. For the volume examined in Fig. C.1, f₆₀ = 0.25. Thus, even relaxing the criterion, accreted stars turn out to dominate only a limited region of the E − L_z space. For this specific volume and this specific simulation, the probability that a given region of the E − L_z space is mostly made of in-situ stars is thus very high, approximately 75%. We have repeated this calculation for all the volumes defined in Fig. 5, for the remnant of the 1 × (1:10) simulation. The corresponding values, f₉₀ and f₆₀ are reported in Table 3. On average, for the volumes located at 8 kpc, f₉₀ = 0.16 and f₆₀ = 0.26. For those located at 12 kpc, the fraction of the E − L_z space dominated by accreted stars increases slightly (f₉₀ = 0.21 and f₆₀ = 0.34), mainly due to the decreasing density of the in-situ halo population.

The analysis presented in the first part of this section concerns the case of a single 1:10 accretion. It is thus important to understand how the results are sensitive to the accretion history of the galaxy and in particular to the number of accreted satellites. For this, we need to examine the cases of 2 × (1:10) and 4 × (1:10) accretions. By examining these two cases, we do not assume, of course, to have explored the full parameter space. But we can, in any case, derive some lessons, as detailed below. Similarly to the analysis above, Figs. C.2 and C.3, respectively, show the distribution in E − L_z space of all halo stars in a solar vicinity volume and that of all in-situ and accreted stars (also separating the contribution of the different satellites) with their relative contributions to the total distribution. As before, halo stars were selected as those satisfying the kinematic condition km s^-1. There are some points worth retaining from this analysis:

• 1.

  in all cases, neither the total distribution of halo stars, nor that ofaccreted stars, nor that of in-situ stars is smooth in the E − L_z space;

• 2.

  the higher the number of accretion events, and thus the accreted stellar mass, the larger the distribution of in-situ stars in E − L_z space;

• 3.

  each satellite redistributes over a large area of the E − L_z space in the solar vicinity, giving rise to several substructures;

• 4.

  stars initially belonging to different satellites, even with significantly different initial orbits, as is the case of the 4 × (1:10) merger, can significantly overlap in the E − L_z space, to the point at which the global distribution of accreted stars is barely decipherable.

As a consequence,

• because of Point 1, substructures in this space cannot beinterpreted as an indication of the extragalactic origin of the starsforming them;

• because of Point 2, it is not possible to define a priori what the extension of in-situ halo stars in the E − L_z space should be: depending on the specific accretion history of the galaxy, in-situ stars heated by the interaction can occupy a variable extent of this space. Any cut at any particular value of E and L_z, as tentatively suggested by Ruchti et al. (2014, 2015), is arbitrary. Moreover, as we discuss more extensively in Sect. 4, the in-situ halo stellar population in our model is only made up of stars initially in the disc and then kinematically heated by the interaction to form the halo. In other words, in our models, an initial population of in-situ halo stars is missing. If this population was taken into account, the distribution of the in-situ population in the E − L_z space would perhaps become even more extended, depending on, among other factors, the intrinsic amount of rotation this initial in-situ halo population would have. Modeling this initial population is beyond the scope of this paper, but it is worth taking into account that throughout this analysis we may probably still underestimate the role of in-situ halo stars;

• because of Point 3, from the number of substructures found in this space we cannot trace back the number of accretion events experienced by the galaxy;

• because of Point 4, it is not possible to separate the different satellite progenitors in this space.

Finally, we can repeat the analysis made for the 1 × (1:10) merger and measure which fraction of the E − L_z space is dominated by accreted stars in the case of our 2 × (1:10) and 4 × (1:10) mergers. The results of this analysis are reported in Tables 4 and 5, respectively. It is only for the 4 × (1:10) simulation that the fraction of the E − L_z space is significant, with an f₆₀ ~ 0.5 and f₉₀ ~ 0.3, on average. These higher fractions are mostly due to the presence of a significant fraction of accreted stars on retrograde orbits due to the accretion of satellite #4 and partially also of satellite #3 (see also Fig. C.3). Even in this case, however, the role of the in-situ population is still important, half of the E − L_z space still being dominated by in-situ stars.

3.2. On the L_⊥ – L_z space

Differently from the energy, the determination of the angular momentum components does not require any knowledge of the Galactic potential. Moreover, L_z is conserved in an axisymmetric potential and L_⊥, even if not strictly conserved, may vary only marginally (see Sect. 2, Chap. 3 in Binney & Tremaine 1987). This is the reason why the L_z − L_⊥ space has been used several times as a natural space to look for the presence of stellar streams, ever since seminal work by Helmi et al. (1999). By studying a sample of approximately 100 metal-poor red giants and RR Lyrae stars within 1 kpc of the Sun, they pointed out that while at values of L_⊥ ≤ 1000 kpc km s^-1 and | L_z | ≤ 1000 kpc km s^-1, the distribution appears relatively symmetric, with stars spanning all values of L_z. For larger L_⊥, the distribution appears asymmetric with respect to L_z, with an excess of stars on prograde (i.e., in the direction of Galactic rotation) orbits, showing up as an over-density. According to their work, indeed, for such high values of L_⊥, very few stars appear on retrograde orbits (L_z< 0) and there is a lack of stars on polar orbits (L_z = 0). By comparing the observed distribution of stars in the L_z − L_⊥ space with Monte Carlo simulations of smooth, non rotating stellar halos, the over-density found at (L_z,L_⊥) ≃ (1000, 2000) kpc km s^-1 was shown to be significant and was interpreted as evidence of an accreted stream, known as the “Helmi stream” after this discovery. Other studies have subsequently confirmed the presence of this overdensity in the L_z − L_⊥ space, confirming its extra-galactic nature (Chiba & Beers 2000; Re Fiorentin et al. 2005; Dettbarn et al. 2007; Kepley et al. 2007; Klement et al. 2009; Smith et al. 2009; Re Fiorentin et al. 2015) and concluding that approximately 5% of the local stellar halo should be made of stars belonging to this stream. Apart from the Helmi stream, possible others have potentially been detected in the L_z − L_⊥ plane (see, for example Kepley et al. 2007; Smith et al. 2009; Re Fiorentin et al. 2015). A point that is particularly critical in the search for streams and substructures in the L_z − L_⊥ space is the comparison that is often made with what is supposed to be the distribution of in-situ halo stars in this space (see, for example Helmi et al. 1999; Kepley et al. 2007; Smith et al. 2009): in-situ halo stars are assumed to have a smooth distribution in the L_z − L_⊥ space and to not rotate. This is a critical assumption, which affects the significance of any possible detection.

We thus concentrate the following analysis on the investigation of the locus occupied by in-situ halo stars in this plane and on understanding and quantifying their overlap with accreted stars. Once again, we remind the reader that the stellar halo, in our simulations, is entirely made of disc stars kinematically heated by the interaction. No modeling of a stellar halo pre-dating the accretion event is included in this work.

3.2.1. Looking for streams in the solar vicinity: accreted, in-situ stars and their overlap

In Fig. 7, we show the distribution of halo stars, for some of the solar vicinity volumes defined in Fig. 5 and for the case of the 1 × (1:10) accretion, the 2 × (1:10), and the 4 × (1:10), respectively. Halo stars were selected kinematically, using the same criterion already adopted in Sect. 3.1.3 and in Figs. C.1–C.3: among all stars in the chosen solar volume, halo stars are those satisfying the condition km s^-1.

The main points we want to make from this analysis are the following:

• 1.

  Halo stars occupy a large portion of the L_z − L_⊥ space, whose extension depends on the number of accretion events experienced by the galaxy, and on the orbital characteristics of the accreted satellites (prograde or retrograde orbits). In our models, it is the 4 × (1:10) accretion that leads halo stars to redistribute over the largest area of the L_z − L_⊥ space.

• 2.

  In general, the distribution is asymmetric with respect to the axis L_z = 0, towards positive L_z values (i.e., prograde motions). The strength of this asymmetry depends on the accreted mass – the higher the number of accretion events, the more symmetric the distribution becomes – and on the presence of satellites on retrograde orbits. The cause of this asymmetry cannot be attributed solely to the presence of accreted stars. Indeed, as Fig. 7 shows, in-situ halo stars also show an asymmetric distribution, skewed towards positive L_z. Again, the strength of this asymmetry depends on the number of accretion events and on their characteristics – also for in-situ stars, the distribution tends to become more symmetric with respect to the L_z = 0 axis in the case of the 4 × (1:10) interaction – but in all cases, at high values of L_⊥, the distribution of in-situ halo stars becomes skewed towards more prograde, inclined orbits. Even if the accretion event causes an important redistribution of the initial angular momentum of the main galaxy and in particular of its z −component, to the point that some in-situ stars end up on retrograde orbits, overall, in-situ halo stars retain part of their initial rotation, forming a slow-rotating stellar halo. This finding is fundamental: it implies that any observational evidence of an asymmetric distribution in the L_z − L_⊥ space skewed towards prograde, inclined orbits is not in itself an indication that this region of the space is dominated by accreted stars. Any accretion event indeed generates an in-situ halo population (made of pre-existing disc stars heated by the interaction), whose distribution, at large values of L_⊥, is skewed towards large values of L_z.

• Overall, the overlap between the in-situ and accreted population in the L_z − L_⊥ plane is considerable everywhere. There is no particular region where the contribution of accreted stars appears dominant. For values of (L_z − L_⊥) approaching (0,0), the in-situ population dominates. While moving towards the periphery of the distribution, the two populations show comparable stellar densities overall. This is visible, for example, in the rightmost panels of Fig. 7, where the density contours of the in-situ and accreted populations are reported: they show a striking similarity everywhere, suggesting that even in the periphery of the distribution (for large values of L_z and L_⊥, for example), the probability of finding in-situ stars is high. To better quantify this point, for each solar vicinity volume, we have selected two regions of the L_z − L_⊥ plane: one corresponding to the region where the Helmi stream has been found observationally (i.e., in units of 100 kpc km s^-1, 5 <L_z< 15 and 14 <L_⊥< 25, hereafter referred to as the Helmi region, see Kepley et al. 2007) and one characterized by similar values of L_⊥ but more extreme values of L_z (15 <L_z< 25, hereafter referred to as the extreme prograde region). For each of these two regions, for each solar vicinity volume and for all the three simulations, we have evaluated the number of satellite stars contained in that region and compared it to the total number of stars found in the same region. These values are reported in Tables 6–8, for the 1 × (1:10), 2 × (1:10), and 4 × (1:10) simulations, respectively. While the fraction of satellite stars depends on the volume under consideration, we found that, on average, there is no clear dependence of this value on the number of accreted satellites and most importantly that in all cases a significant fraction of in-situ stars is present in the selected regions, even when the Milky Way-type galaxy experiences multiple accretion events. Averaging over all the volumes, in the Helmi region, the average fraction of accreted stars is 0.39 in the case of the 1 × (1:10) simulation, 0.37 for the 2 × (1:10) simulation, and 0.5 for the 4 × (1:10) simulation. For the extreme prograde region, these fractions are 0.63, 0.46 and 0.22, respectively.

• 4.

  Finally, in reference to the smoothness of the L_z − L_⊥ space, the overall distribution does not appear to be smooth, thus confirming observational findings (see, for example Helmi et al. 1999; Smith et al. 2009; Re Fiorentin et al. 2015). The presence of substructures is visible among accreted halo stars, in agreement with previous models (Helmi et al. 1999; Kepley et al. 2007; Re Fiorentin et al. 2015). Most importantly, however, a result, to our knowledge, thus far undiscovered, the in-situ halo population is also non-uniformly redistributed in the L_z − L_⊥ space. Some examples of this clumpy, in-situ halo distribution can be appreciated in Fig. 7. Over-densities of in-situ stars appear not only at low values of L_⊥, thus for low inclination orbits, but also at values as large as L_⊥ ~ 20. This suggests that extreme caution must be taken in interpreting the origin of substructures and clumps in the L_z − L_⊥ space: the presence of isolated/dense groups of stars in some regions of this space, or of groups with extreme prograde or retrograde orbits, is not in itself a probe for extragalactic origin of these stars.

3.2.2. In-situ and accreted stars in velocity spaces

Before discussing the signature of accretion events in the APL space, it is worth commenting on the distribution of in-situ and accreted stars in velocity spaces.

The velocity of stars belonging to the Helmi stream was investigated in Helmi et al. (1999), Kepley et al. (2007) (see also Klement et al. 2008; Helmi 2008; Smith 2016). In the v_r − v_z (i.e., Galactocentric radial velocity – vertical velocity) plane, stars in the stream appear as two separate groups, with v_z clumping at approximately ±200 km s^-1 and radial velocities distributed in the interval − 200 km s^-1 ≲ v_r ≲ 200 km s^-1. The split found in v_z and visible also in the v_z − v_φ plane, where v_φ is the tangential velocity, can be understood with the presence of streams of stars having similar L_⊥, but moving in opposite directions (i.e., inward and upward) with respect to the Galactic plane. These observations were compared to N-body simulations aimed at modeling the accretion of a satellite galaxy in an analytic Milky Way potential in Helmi et al. (1999), Kepley et al. (2007), Helmi (2008). It was found that the properties of stars in the Helmi stream could be reproduced if the stream was the remnant of a satellite with an initial internal velocity dispersion of 12 − 18 km s^-1, probably resembling the Fornax or the Sagittarius dwarf galaxy (Helmi et al. 1999) and accreted between 6 and 9 Gyr ago (Kepley et al. 2007).

The message one can retain from these works is that information can be retrieved by looking solely at the kinematics of stars at the solar vicinity: their in-situ or extragalactic origin and, in the latter case, the mass of the progenitor satellite which generated the stream. Our models delineate a more complex scenario.

The v_r,v_φ, and v_z velocities of simulated halo stars in some of the solar vicinity volumes previously studied are shown in Figs. 8–C.5. As in Sects. 3.1.3 and 3.2, for each solar vicinity volume in Fig. 5, we have selected halo stars as those satisfying the condition km s^-1, and, among those, stars in the Helmi region are defined as those with 5 <L_z< 15 and 14 <L_⊥< 25 (in units of 100 kpc km s^-1, see previous section and the rightmost panels of Fig. 7). Their velocities are reported in Figs. 8–C.5 for some of the examined solar vicinity volumes and for all the simulations. The velocities of all simulated stars occupying the Helmi region are reported in the top panels of each of these figures and compared to the velocities of all stellar particles in the selected volume. One can see that the kinematic characteristics found for the Helmi stream at the solar vicinity (Helmi et al. 1999) are recovered: our modeled Helmi stars show extreme values of v_z and v_φ when compared to the whole set of halo stars at the solar vicinity, they show the split in the v_z distribution also found in the observations (see panels in the first and last columns of Figs. 8–C.5) and their v_r and v_φ values are similar to those attained by the observed sample. However, the finding that the modeled stars reproduce the kinematic distribution of observed ones does not say anything about the nature of these stars (in-situ or accreted). The central and bottom panels of Figs. 8–C.5 indeed show that in-situ and accreted stars in the Helmi region have the same velocities: similar v_r, v_φ, and v_z. In particular, in-situ stars also show the characteristic split in v_z observed for stars in the Helmi stream.

That in-situ and accreted stars in the Helmi region have the same velocities is not surprising. They belong to the same region of the L_z − L_⊥ plane and as a consequence have similar angular momenta. Further, they have been selected in a limited volume around the Sun, which limits the range of distances they span. The angular momenta and the distances being fixed and similar, this naturally also leads to similar velocities. The conclusion that we draw from this analysis is that all stars in the Helmi region should share the same kinematics, independently from their in-situ our accreted origin. The fact that, under certain conditions, N-body models can reproduce the velocities observed in the Helmi region by means of an accreted stream is not a probe in itself for accreted stars observed in that region: once the spatial volume is fixed, all stars in that volume with similar angular momenta will have similar velocities, independently from their accreted or in-situ nature.

Not only do our models suggest that the use of these velocity spaces cannot solve the question of the in-situ or accreted origin of the stars that make them, but also that it is not possible to recover the mass of the progenitor satellite, if any, from them. While Helmi et al. (1999) and Kepley et al. (2007) suggest that the progenitor of the Helmi stream may have been a galaxy similar to Fornax or Sagittarius, our models show that the kinematic properties observed for this stream can be reproduced with one (or multiple) satellite(s) significantly more massive, each of them having an initial stellar mass approximately 100 times larger than the current stellar mass in Fornax (de Boer et al. 2012).

3.3. On the R_apo–R_peri space

We end our quest for accreted stars by exploring the APL. The distribution of solar vicinity stars from the Nordström catalogue in this space has already been explored by Helmi et al. (2006), who compared the observed distribution with predictions of N-body models of satellites accreted in a fixed analytic Milky Way potential and with cosmological simulations. Since, in both cases, the distribution of in-situ stars, and their possible overlap with accreted populations, has not been investigated, we discuss the issue in this section. In the following, the apocenters and pericenters, indicated as R_apo and R_peri respectively, are the true apocenter and pericenter distances spanned by the particles in the last Gyr of evolution, that is, no assumption on the gravitational potential is made a posteriori to reconstruct the orbital parameters.

Figures C.6 and C.7 show the APL space for the 1 × (1:10) and 2 × (1:10) simulations, respectively, for all stellar particles in a given solar vicinity volume. From these plots we can deduce the following results:

• as for the E − L_z and L_⊥ − L_z spaces discussed in previous sections, the overlap of in-situ and accreted stars is also substantial in the R_apo − R_peri space, to the point that no clear distinction can be made among in-situ and accreted stars on the basis of this space only. Similar conclusions are found for the R_apo − L_z space (not shown).

• for some of the volumes explored, some clumps of stars are found at large apocenter values (R_apo ~ 40 kpc), but the origin of these clumps is not always extragalactic. For example, while in Fig. C.6 (top row), the over-density found at (R_peri,R_apo) = (3, 40) kpc is made of accreted stars, in Fig. C.7 the two evident over-densities at (R_peri,R_apo) = (6, 45) kpc and (R_peri,R_apo) = (8.5, 38) kpc are mostly made up of in-situ stars.

• The global distribution of all the stars in the volume in the apocenter-pericenter space reveals the presence of several streaks, already noticed by Helmi et al. (2006) in their models. Each of these streaks spans a large range of eccentricities, e, from e ≤ 0.1 to e ≥ 0.7. They are visible both in the accreted and in-situ components, both on a large scale (top row panels of Figs. C.6 and C.7) and on a zoomed-in region of the R_apo − R_peri space, defined by 0 ≤ R_peri ≤ 8 kpc and 7 ≤ R_apo ≤ 15 kpc, whose area is similar to that studied by Helmi et al. (2006). As we have verified, these streaks are made up of stars with similar energies, but which span a large area of the vertical component of the angular momentum L_z (see Figs. C.6, C.7). They are associated with spiral-like features, or rings of different sizes, induced by the merger event, similarly to the structures discussed in Gómez et al. (2012). These streaks are not only visible in the stellar (thin and thick) discs, as already noticed by Gómez et al. (2012), but also among halo stars. We highlight those in the bottom rows of Figs. C.6 and C.7 that show the highest velocities ( km s^-1) in the LSR reference frame.

• In the zoomed-in region, according to our models, in-situ stars occupy most of the R_apo − R_peri space, independently from the number of accreted satellites and the solar vicinity volume under study. This is also true for the region lying between the lines of constant eccentricities e = 0.3 and e = 0.5 (see central rows in Figs. C.6 and C.7), which corresponds to the region where over-densities found in the solar vicinity data from the Nordström catalogue have been associated with accreted streams (Helmi et al. 2006). Our models suggest that, in this region, the contribution of in-situ stars heated by the interaction and going on to populate the thick disc-inner stellar halo is substantial.

• Finally, for the over-densities found in the R_apo − R_peri distribution of stars in the Nordström catalogue that follow a diagonal pattern at eccentricities 0.3 ≤ e ≤ 0.5 and for which an extragalactic origin has been suggested, a trend has been reported by Helmi et al. (2006): these stars tend to have progressively larger L_z as the R_apo and R_peri increase. We notice from the bottom panels of Figs. C.6 and C.7 that, in our models, this trend (an increase of L_z with R_apo and R_peri ) is found among both accreted stars and the in-situ population.

To conclude, both the diagonal distribution of stars in the R_apo − R_peri space and the presence of streaks or over-densities and the trend observed for the angular momentum L_z with R_apo and R_peri, are features common both to accreted and in-situ stars. Without detailed chemical abundances, the differentiation of the extragalactic or in-situ origin of stars in this space is, according to these results, unfeasible.

Similar results are found also for the 4 × (1:10) simulation.

3.4. Leaving the solar volume: kinematic detection of streams on a 10 kpc scale

Most of the previous analysis, except that presented in Figs. 2–4, has been performed selecting stellar particles in restricted solar volumes (3 kpc in radius). These are the typical volumes over which the kinematic search of stellar streams has been focused until now (Helmi et al. 1999; Kepley et al. 2007; Morrison et al. 2009; Smith et al. 2009; Re Fiorentin et al. 2005, 2015). In the coming years, with Gaia and spectroscopy-related surveys such GALAH, WEAVE, and 4MOST, we will have an almost complete view of stars up to ~3 kpc from the Sun: distances, kinematics, ages and detailed abundances. This full set of quantities should thus provide enough information to separate patterns and reconstruct the origin and the formation sites of “local” stars. However, beyond the ~3 kpc-sphere, at larger distances from the Sun, uncertainties in dating stars will possibly be too large to determine ages with sufficient accuracy² and detailed abundances from large spectroscopic surveys with a precision of ~0.03–0.05 dex will also be difficult to obtain³. It is thus at those scales (typically between 3 and 10 kpc from the Sun), that it is fundamental to understand if the search for the fossil records of the Galaxy by kinematics alone is feasible and meaningful.

In Fig. C.8, the analysis of these extended solar volumes, up to distances of 10 kpc from the Sun, is shown, for the E − L_z plane. For each simulation, we have selected all stars with distances inside 10 kpc from the Sun, where four different Sun positions have been chosen, corresponding to the centers of volumes 1 and 5 in Fig. 5, at distances of 8 kpc and 12 kpc from the galaxy center. We have also excluded from this analysis stars at vertical distances z from the plane below 3 kpc to maximize the fraction of halo stars in the examined samples. In Fig. C.8, only one of these volumes is shown, the conclusion of the analysis being the same for all four volumes examined. This figure shows that all the problems that affect the local samples are also present in a 10 kpc–extended solar volume: significant overlap of in-situ and accreted stars, in-situ stars that dominate a large part of the E − L_z space, and non-smooth distribution for both populations. Making a blind test; looking at any of the distributions shown in the top row of Fig. C.8, it appears impossible to establish how many different satellites have contributed to determining those distributions, the masses of the progenitor systems, and which fraction of the lumpy regions is made up of in-situ and accreted stars, respectively. Similar conclusions are reached for the L_⊥ − L_z space (see Fig. C.9).

4. Discussion

4.1. In-situ halo stars, the elephant in the room

All the results presented in the previous sections strongly suggest a problem: the kinematic search for streams of accreted satellites in the Milky Way cannot be achieved neglecting the presence of the in-situ stellar population. In the simulations presented in this work, as recalled several times in this manuscript, there is no in-situ stellar halo before the interaction(s). The in-situ halo found at the end of the simulations is only the result of the heating of a pre-existing stellar disc and no other channel for the formation of the in-situ halo, such as those described by Cooper et al. (2015), has been taken into account. That heating from satellite accretions is effective in forming or contributing to forming a stellar halo has already been eluded to in a number of publications Zolotov et al. (2009), Purcell et al. (2010), Font et al. (2011), Qu et al. (2011), McCarthy et al. (2012), Cooper et al. (2015). Here, however, we show that the consequence of this heating is important for any search in integrals-of-motion or kinematic spaces:

• Heated halo stars are not smoothly distributed: this result is incontradiction with the usual assumption made in the literaturethat the in-situ halo should be smooth in those spaces. At least thepart of the in-situ halo that results from satellite heating is, on thecontrary, structured, even several Gyr after the accretion hastaken place. Halo stars with a clustered distribution in kinematicor integrals-of-motion spaces are not necessarily accreted. Theprobability that these clumps lie in regions dominated by in-situstars is indeed very high (see Tables 3–8).

• Heated halo stars rotate: this result is fundamental for all studies that look for accreted streams in the L_z − L_⊥ space and discriminate between the in-situ and the accreted populations by assuming that the in-situ halo should not rotate. Our models indeed show that this assumption, although controversial from the observational point of view (see, for example, Chiba & Beers 2000; Carollo et al. 2007; Kepley et al. 2007; Smith et al. 2009; An et al. 2015), is not valid at least for the part of in-situ halos made up of heated thin/thick disc stars. As we show in Fig. 7, indeed, the distribution of in-situ halo stars in the L_z − L_⊥ space presents an excess of prograde orbits (i.e., positive L_z). As a consequence, the distribution of those in-situ halo stars is not symmetric with respect to the L_z = 0 axis. This asymmetry has been found in observation samples of stars at a few kpc from the Sun for prograde orbits with high values of L_⊥. Stars in these regions have been interpreted as accreted stars precisely because it was assumed that any in-situ population should be non-rotating and thus symmetric with respect to the L_z = 0 axis. Contrary to this common assumption, here we show that in-situ heated halo stars are expected to rotate. The amount of rotation depends on the number of accretion events and on the total accreted mass: the larger the accreted mass, the slower the in-situ halo rotates.

• Heated stars overlap with accreted stars: this finding affects all the spaces studied in this paper; E − L_z, L_⊥ − L_z, and APL, and the proportion of in-situ stars is so important in each of these spaces (see, for example, Tables 3–8 and Figs. 4, 6–C.7) that it leads us to conclude that the search for accreted streams in those spaces is mostly inefficient, simply because the probability of finding in-situ stars in a given region of those spaces is significant everywhere and it is not possible to define a priori where the chance of finding accreted stars is the highest. The extension and location of the region dominated by accreted stars depend indeed on a number of unknown parameters; parameters that we would like to constrain by a search in kinematic spaces, such as the number of accreted satellites, their masses and their orbital properties.

To the difficulty in separating in-situ and accreted stars, we need to add the difficulty in recovering the properties of the accreted satellites (in particular their number and their mass) from their clumpy distribution in kinematic and integral-of-motion spaces. For example, Fig. 3 shows that, even in the ideal case where we are able to separate in-situ from accreted stars, the number of ripples and clumps found for the accreted population in the E − L_z space does not reveal anything about the number of progenitors beyond this distribution: a simple 1:10 merger is able to produce tens of distinct structures in this space, each of these with sizes and densities similar to those generated by a 1:100 accretion (see below). This finding is common to all the analyzed spaces. Another example is given by the v_r − v_φ, v_z − v_φ, and v_r − v_z spaces: on the basis of the velocities found for accreted stars in the solar vicinity (their values and characteristics, see Figs. 8, C.4, and C.5), it is not possible to derive the number of accretion events which contributed to shaping the accreted population.

4.2. Chemistry and ages, the key for differentiating the origin of stars in the Galaxy

In attempting to overcome this degeneracy, one may be tempted to think that a simple cut in metallicity (for example, selecting only stars with [Fe/H] ≲ − 1 dex) would be sufficient to eliminate the contamination of in-situ stars. The reasoning behind this choice would be to think that heated halo stars previously in the disc should have metallicities typical of the Milky Way disc, that is [ Fe / H ] ≥ − 1 dex. Following this reasoning, a cut in metallicity that removes the metallicities of current disc stars would allow for removal of in-situ contamination. However, this reasoning risks being incorrect because the metallicities expected for the heated halo stars would be the metallicities of the disc at the time the accretion events took place and not the metallicities of the current (i.e., redshift z = 0) Milky Way disc. Thus, the metallicities of the in-situ heated halo population critically depend on the exact moment when relatively massive mergers (mass ratio ≲1:10) took place. In the cosmological simulations analyzed by Tissera et al. (2012), for example, in-situ and accreted inner halo stars have similar metallicities (see Table 5, in their paper), generally below [ Fe / H ] = − 1 dex. Any cut on the basis of metallicity only would thus not be sufficient to remove the in-situ population, simply because part of this population has possibly been kinematically heated in the halo at early times when the most massive mergers took place and thus has metallicities typical of current halo stars.

To solve the question of the origin, site, and mode of formation of stars in the Galaxy, we need to look for detailed chemical abundances and/or ages. ΛCDM models indeed predict that several chemical and chemo/age patterns should be expected in Milky Way-like galaxies, for stars originating in different environments (see Gibson et al. 2013; Few et al. 2014; Snaith et al. 2016, for some recent works). However, the ability to find and differentiate those patterns in future spectroscopic samples of stars crucially depends on the precision of the abundance estimates achieved by future surveys (see Figs. 12 and 13 in Snaith et al. 2016). From a solar vicinity sample of stars, where abundances with precisions below 0.1 dex have been attained, it has been possible to reveal a number of chemical (or chemical/age) patterns. The dominant patterns seem characteristic of the evolution of the in-situ galactic stellar populations (Haywood et al. 2013; Spina et al. 2016; Nissen 2016), with possible dilution from gas accretion in the outer disc (Haywood et al. 2013; Snaith et al. 2015), while some others possibly have an extragalactic origin (Nissen & Schuster 2010). This approach seems very promising and its predictive power, in terms of ability to separate populations of different origins, seems to surpass that of kinematic diagnostics, as can already be evinced from stellar samples at the solar vicinity (see below).

4.3. An illustrative example: the mapping of a chemically selected, local stream into kinematic spaces

To illustrate the potentiality of the search for accreted/in-situ populations when detailed chemical abundances are available, in Fig. 9 we show the distribution of halos and thick disc stars in the solar vicinity in the E − L_z, L_⊥ − L_z and R_apo − R_peri spaces from the sample studied by Nissen & Schuster (2010) and subsequently in Nissen & Schuster (2011, 2012), Schuster et al. (2012). We chose this sample because by studying detailed chemical abundances, Nissen & Schuster (2010) show that the halo stars that form it can be separated into two distinct populations on the basis of their [α/Fe] content (see Fig. 1, Nissen & Schuster 2010). The high-α sequence appears to be in chemical continuity (with some overlap) with the thick disc and, as proposed by Nissen & Schuster (2010), it may be made of stars originally in the disc or bulge of the Galaxy and subsequently heated by merging satellites. The currently favored hypothesis for the low-α sequence is that it is made of stars accreted from dwarf galaxies, with some of them possibly associated to ω Cen progenitor. This sample thus very probably contains both an in-situ and an accreted halo population and since parallaxes, radial velocities, and proper motions are available (see Nissen & Schuster 2010; Schuster et al. 2012, for details), we can derive orbital properties and make use of them as a natural testbed for probing the efficiency of the detection of accreted streams by kinematic and integrals-of-motion diagnostics.

To achieve this aim, we have integrated the orbits of all stars in this sample in the Allen & Santillan (1991) Galactic potential for 7 Gyr. We show the value of their energies, angular momenta, R_apo, and R_peri in Fig. 9. This figure clearly shows that in none of those kinematic spaces does the sequence of low-α accreted halo stars show up as a cluster or substructure distinct and separated from the whole sample. As already demonstrated by Schuster et al. (2012), it is true that the low-α sequence shows, on average, larger values of energies and R_apo than those reached by the high-α, in-situ sequence, but while the kinematics can be used to study differences among those (chemically identified) populations, it is clear that it cannot be used to identify the accreted stream and separate it from the bulk of the in-situ population. It would not be possible to recognize and separate the two halo populations in any of the spaces shown in Fig. 9 without marking the different populations with different colors. There is a non-negligible overlap between the accreted and in-situ populations: this is particularly striking in the L_⊥ − L_z plane, where all the accreted stars overlap with the in-situ population, except for two stars (i.e., G112-43 and G112-44, a couple of visual binaries. See Nissen & Schuster 2010; Allen et al. 2000, for their chemical properties and orbit), which lie in the upper-right part of the diagram, at values of L_z and L_⊥ comparable to those of the Helmi region.

This example thus illustrates the (generally) inefficient power of kinematic diagnostics in finding accreted streams and, at the same time, how promising the search for accreted streams in chemical abundance spaces can be when compared to kinematic detections (cf. Fig. 1 in Nissen & Schuster 2010, with Fig. 9 in this paper). The power of this approach is considerable and is even expected to increase once ages and detailed chemical abundances of several hundreds thousand stars are available. However note that for this, we need very high quality spectroscopic data, the separation in [α/ Fe] between the two halo sequences found by Nissen & Schuster (2010) being typically less than 0.1 dex.

4.4. The kinematic signature of mergers with a mass ratio ~1:100

Mergers with mass ratios similar to those analyzed in Sect. 3 appear inevitable for Milky Way-type galaxies in ΛCDM simulations (Read et al. 2008; Stewart et al. 2008; De Lucia & Helmi 2008; Deason et al. 2016). However, these simulations also predict that accretions of satellites with smaller masses should also take place. Here we want to understand whether it would be possible to recover these low-mass accretions in kinematic or integrals-of-motion spaces more easily than what our models suggest for more massive satellites. In the E − L_z space, for example, we have seen that the non-conservation of energy and L_z is mostly driven by dynamical friction. Because the latter depends on the mass of the accreted galaxy, one would expect low-mass satellites to conserve their initial location in the E − L_z space, without experiencing the significant changes in energy and angular momentum quantified for the more massive systems in Sect. 3.1. Motivated by this question and by the need to quantify the ability to recover accretions with low-mass ratios in a Milky Way-type galaxy by means of kinematics alone, we ran some supplementary simulations, where the Milky Way-type galaxy, after experiencing a 1 × (1:10) or a 2 × (1:10) merger, accretes four or eight satellites with mass ratio 1:100. More precisely, we took the final snapshot of the simulations 1 × (1:10) and 2 × (1:10) analyzed in Sect. 3 and for each of these two new “initial conditions” we ran two simulations, accreting, 4 or 8 satellites, respectively, with masses one hundred times smaller than those the Milky Way-type galaxy initially had before the 1:10 mergers. The low-mass satellites are re-scaled versions of the Milky Way-type galaxy: their masses are a factor 100 smaller, and their sizes have been divided by a factor of 10. For computing-time reasons, for each of these simulations, we reduced the number of particles of both the Milky Way-type galaxy and its satellites by a factor of 10 and re-scaled the masses of the single particles accordingly. This means, for example, that the 1 × (10) + 4 × (1:100) merger contains approximately 2.8 million particles in total, instead of 28 million particles, expected if we had used the resolution adopted for all simulations described in Sect. 3. Each low-mass satellite is made of 25 000 particles and does not contain any population of globular clusters. The lower resolution employed to run these simulations does not allow for appreciation and study of the fine structure observed in the integrals-of-motions spaces discussed in Sect. 3. However, it is still possible to quantify the contamination and overlap between in-situ stars and the different accreted stellar populations. An example of the difficulty in deciphering the signature left by low-mass mergers is given in Fig. C.10. In this example, the Milky Way-type galaxy, after experiencing two mergers with a mass ratio 1:10, accretes eight satellites with mass ratios 1:100. The distribution of all stars in the E − L_z space is shown here, independently of their location in the galaxy. The points that we want to emphasize from the analysis of these supplementary simulations are the following:

• as found in the previous section, even in the case of a2 × (1:10)+8 × (1:100) merger, such as that shown in Fig. C.10, most of the E − L_z space is dominated by in-situ stars. The only region where the contribution of satellites becomes dominant is that at energies E ≳ − 2. Note, however, that this region (dominated by accreted stars, i.e., E ≳ − 2) would be inaccessible if we had restricted the search to a sphere of ~10 kpc from the Sun.

• Low-mass satellites are characterized by structures, in the E − L_z plane, comparable to one of the several generated by a more massive accretion. Indeed, if we compare the over-densities generated by 1:100 mass-ratio mergers, they are similar in sizes and densities to one of tens of clumps generated by a 1:10 merger (cf., for example the distribution of the satellites 1 and 2 with mass ratio 1/100 to part of the distribution generated by a 1:10 merger).

To conclude, even if, in general, mergers with mass-ratios of approximately 1:100 tend to better preserve their initial distribution in the E − L_z space, for the Milky Way to experience even one 1:10 merger in its past would be enough to render the E − L_z space barely decipherable. Mergers with such mass ratios indeed erase the signature of the less massive ones, because they dominate in mass and because each of the over-densities they generate in the E − L_z space is indistinguishable from those generated from less massive accretions.

4.5. Limitations of the modeling presented in this work

Before moving to the main conclusions of this work, it is necessary to recall its main limitations, which we wish to address in future studies, and to discuss some possible avenues for overcoming them.

One limitation comes from the limited range of orbital parameters and internal parameters (such as: satellite to MW-type-galaxy-mass ratios, M/L ratios and relative number of particles for the baryonic and dark matter components, etc.) explored. Even if a study with more extended parameters space is suitable, with initial conditions possibly taken from cosmological simulations, we think the current choice does not impact the main results of this work, that is, the inner halo formation by disc heating and the overlap of in-situ and accreted stars in kinematic and integrals-of-motion spaces. Even if not specifically focused on integrals-of-motion spaces, several works have investigated the impact of satellite accretions on heating pre-existing stellar discs and the resulting kinematic properties of in-situ and accreted populations (a non-exhaustive list is given in Sect. 3.1.2). The results presented in those papers have explored a large interval of internal and orbital parameters, showing in particular that the fraction of satellite stars over the total number of stars found at a given height above the plane of the remnant galaxy depends very weakly on the orbital conditions of the merger (see Fig. 16 in Villalobos & Helmi 2008). This gives us confidence that even if we have simulated only a few orbits, our results do not depend strongly on the choice of the orbital parameters. We also note that the thickness of the final disc is mostly non-dependent upon the prograde/retrograde orbit (again, see Villalobos & Helmi 2008). On the basis of these results, we feel confident that this point does not constitute a major limitation of our conclusions.

A second limitation of our modeling consists of the absence of an in-situ halo population in the Milky Way-type galaxy that pre-exists the accretion event(s). Indeed, as we have mentioned several times in this paper, in our simulations, the in-situ halo is exclusively made of in-situ disc stars heated by the interaction, and no in-situ halo component pre-exists the accretion(s). It will be thus important to introduce one in future simulations, to quantify its response to the accretion(s), its distribution in integrals-of-motion spaces and the degree of substructures it can show in those spaces. The addition of this component to the modeling will naturally require exploration of a number of different parameters (flattening, rotation, and concentration, among others), an understanding of their impact on the final result.

A third limitation, which is probably the most urgently in need of investigation, is the lack of a gaseous component both in the Milky Way disc and in the satellite(s). It is not easy to anticipate how the inclusion of gas in the simulations can impact the results. It has been shown (Moster et al. 2010; Qu et al. 2011) that the presence of a thin gaseous disc can limit the amount of disc heating: for gas fractions of the order of 20% (40%), the heating may be reduced by 25% (40%). For accretions occurred in the early history of the Galaxy, typically before redshift z ~ 1, the fractions of gas in the Milky Way disc may have been significantly higher (with f_gas ≥ 50%) than those simulated in the above cited works. Moreover, in those early phases of the Galactic evolution it is not even clear how thin the gaseous disc may have been. While in nearby disc galaxies gaseous discs indeed have scale heights between a few tens and 100 pc and velocity dispersions of the order of ~10 km s^-1. At higher redshift gaseous discs of approximately 1 kpc (corresponding to velocity dispersions of the order of 100 km s^-1) have been observed in the continuum and in the gas component (Elmegreen & Elmegreen 2006; Epinat et al. 2012). In the N-body models cited above, the typical scale heights employed for the gaseous component are between 200 and 400 pc. It is thus crucial to quantify the heating of stellar discs and the efficiency in forming an in-situ halo also in more extreme conditions than those modeled so far, also because galactic stellar archeology seems to suggest that the Milky Way may have experienced an intense phase of star formation and significant gas turbulence at those epochs (Haywood et al. 2013; Snaith et al. 2014; Lehnert et al. 2014). If the inclusion of a gaseous component in the Milky Way disc seems inevitable for the simulation of accretion episodes that occurred at early times, the same is valid for modeling the accreted satellites. In the case of gas-rich satellites, the stellar fraction of the baryonic mass would be lower than 100%, as assumed in this paper and a certain fraction of the gaseous mass may be lost by tidal effects/ram pressure stripping before being accreted and converted into new stars. This would naturally lead the accretions to be less “stellar-rich” than those simulated in the present paper, thus effectively reducing the amount of stars of extra-galactic origin that a satellite of a given mass can bring into the inner halo (i.e., inside 20 kpc) of a Milky Way-type galaxy. Finally, it is worth emphasizing that we are using idealized simulations of local mergers to quantify the response of a stellar disc to merger events that may have occurred several Gyr ago, when both the stellar masses of the Milky Way and of its satellite(s) were significantly different from those employed here. Because the properties of the Milky Way at higher z are still largely unknown (in terms of mass and effective radius, etc.) and similarly, those of the satellite galaxies, idealized simulations of local mergers are currently still one of the only viable ways to quantify the response of a stellar disc to accretion events of a given mass ratio, even for those occurring at higher z (with the limitations previously discussed, such as the lack of a gaseous component). Some N-body models similar to those analyzed in this paper have tried to mimic mergers at higher redshift, finding that the evolution and thickening of the disc, and its orbital parameters, do not depend on the redshift, but only on the mass-ratio of the merger (cf., for example, Figs 10 and 14 in Villalobos & Helmi 2008). This is robust evidence of the fact that what matters is the relative mass ratio of the interacting galaxies (see also Quinn et al. 1993; Bournaud et al. 2005, 2009; Qu et al. 2010), and the relative variation of their parameters, rather than the absolute values. As a consequence, at first approximation, local mergers can be re-scaled to also mimic mergers at higher z. This, we suggest caution to be taken due to the limitations of these models, such as, for example, the lack of gas, that can play an important role in the process.

5. Conclusions

In this paper, we have analyzed high-resolution, dissipationless simulations of a Milky Way-type galaxy accreting one or several (up to four) satellites with mass ratio 1:10. These simulations have been complemented by four simulations, at lower spatial resolution, of a Milky Way-type galaxy, accreting several (from four to eight) satellites with mass ratios 1:100, after experiencing one or two mergers with mass ratios 1:10. The novelty of this work with respect to the majority of those already available in the literature, comes from our analysis of fully consistent models, where both the satellite(s) and the Milky Way galaxy are “live” systems, which can react to the interaction, experience kinematical heating, tidal effects and dynamical friction. We have analyzed this set of simulations to investigate the possibility of making use of kinematics information only to find accreted stars in the Galaxy, remnants of past accretion events, that have since then lost their spatial coherence. In particular, we have investigated integrals-of-motion spaces, such as the E − L_z and the L_⊥ − L_z spaces and kinematic spaces such as the R_apo − R_peri and velocity spaces, to understand if we can realistically make use of those spaces to search for accreted streams and if this search is really efficient and meaningful. Our main conclusions are the following.

In the E − L_z space:

• because the energy and angular momentum of a satellite are notconserved quantities during an interaction, each satelliteprovides an origin to several independent over-densities. This isparticularly true for satellites with mass ratios of approximately1:10, which can experience severe dynamical friction during amerger event. In some cases, satellites with a mass ratio 1:100 canalso give rise to several, independent over-densities;

• multiple satellites overlap;

• in-situ stars affected by the interaction(s) tend to progressively populate a region of lower angular momentum and higher energy than the one initially populated (i.e., before the interaction);

• most of the accreted stars overlap with in-situ stars. This point is valid even if the search is restricted to halo stars only, because in our simulations, a substantial part of the inner stellar halo (distances from the galaxy center less than 10–20 kpc) is made of in-situ stars, originally in the disc and then heated by the interaction(s).

In the L_⊥ − L_z space:

• the overlap between the in-situ and accreted population isconsiderable everywhere. There is no particular region where thecontribution of accreted stars appears dominant;

• the distribution is asymmetric with respect to the axis L_z = 0 and skewed towards positive L_z values (i.e., prograde motions). In particular, in-situ halo stars, kinematically heated by the interaction(s), preserve part of their initial rotation. In the L_⊥ − L_z plane, this rotation reveals itself as an excess of stars on prograde orbits. The distribution of in-situ halo stars is thus asymmetric in the L_⊥ − L_z space, similarly to the distribution of accreted stars. This finding is fundamental: it implies that any observational evidence of an asymmetric distribution in the L_⊥ − L_z space skewed towards prograde, inclined orbits is not in itself an indication that this region of the space is dominated by accreted stars. Any accretion event indeed generates an in-situ halo population (made of pre-existing disc stars heated by the interaction), whose distribution, at large values of L_⊥, is skewed towards large values of L_z.

In velocity spaces:

• for a selected spatial volume around the Sun and for similar valuesof L_z and L_⊥, it is not possible to differentiate in-situ from accreted stars;

• in particular, the kinematic characteristics observed for stars in the solar vicinity in the Helmi region can be reproduced by accreted stars and by in-situ stars as well. The finding that, under certain conditions, N-body models can reproduce the velocities observed in the Helmi region by means of an accreted stream does not a clear indication that the stars observed in that region are accreted: once the spatial volume is fixed, all stars in that volume with similar angular momenta will have similar velocities, no matter what their nature, accreted or in-situ;

• not only do our models suggest that the use of these velocity spaces cannot solve the question of the in-situ or accreted origin of the stars that make them, but they also suggest that it is not possible to recover the mass of the progenitor satellite, if any, from them. Our models show that the kinematic properties observed for the Helmi stream can be reproduced with one or multiple satellite(s) significantly more massive than that suggested by Helmi et al. (1999), Kepley et al. (2007) as the progenitor of the stream. The problem of reconstructing the orbital and internal properties of the accreted satellite(s) in the Helmi region by kinematics alone is thus degenerate.

In the APL space:

• as for all the other spaces investigated in this paper, the over-lap of in-situ and accreted stars is also substantial in theR_apo − R_peri space, to the point that no clear distinction can be made between in-situ and accreted stars on the basis of this space only.

• The global distribution of all the stars in APL space reveals the presence of several streaks, already noticed by Helmi et al. (2006) in their models. Each of these streaks spans a large range of eccentricities, from e ≤ 0.1 to e ≥ 0.7. They are visible both in the accreted and in-situ component and, among in-situ stars, not only in the thin and thick discs, as already pointed out by Gómez et al. (2012), but also in the halo.

• all the main characteristics observed among solar vicinity stars in the R_apo − R_peri space (their diagonal distribution in the R_apo − R_peri space, the presence of streakes or over-densities, and the trend observed for the angular momentum L_z, with L_z increasing with R_apo and R_peri ) are common to both the accreted and to the in-situ populations.

In agreement with previous works, we find that all these spaces are rich in substructures, but that the origin of these substructures cannot be determined with kinematics. The in-situ stellar halo, formed as a result of the interaction, is neither smooth nor non-rotating. As a consequence, extreme caution must be employed before interpreting over-densities in any of those spaces as evidence of relics of accreted satellites (see, for example Helmi et al. 1999, 2006; Gould 2003; Kepley et al. 2007; Morrison et al. 2009; Re Fiorentin et al. 2015).

We consider this work to constitute a first step towards a more realistic model of accreted–in-situ populations in kinematic spaces in view of Gaia and its first data releases. It is also intended to highlight the need for caution to be taken during interpretation of current kinematic data where substructures have been detected and an extra-galactic origin is currently favored: all spaces and regions where stellar streams have been detected by means of kinematics are indeed, accordingly to our models, severely dominated by in-situ stars.

On the basis of these simulations, our conclusion is that the kinematic detection of stellar streams in the Galaxy, that is, streams made of stars accreted long enough ago to have lost their spatial coherence, is mostly inefficient. Detailed chemical abundances and/or ages are critical and will definitely be necessary to identify accreted populations, that should appear in those spaces as patterns distinct from those described for in-situ stars.

---

Acknowledgments

This work has been supported by the Ile-de-France Region and the DIM-ACAV, through the grant “Reconstructing the accretion history of the Milky Way through its globular clusters system” and by the ANR (Agence Nationale de la Recherche) through the MOD4Gaia project (ANR-15-CE31-0007, P.I.: P. Di Matteo). The authors also thank the CNRS, for its financial support throughout the MASTODONS project “The origin and evolution of our Galaxy: data validation” (P.I.: F. Arenou). M.M. acknowledges support from the CNR short-term mobility (STM) programme, 2015. We warmly thank: P. E. Nissen and W. J. Schuster, for providing us with a tabular form of the data used in Sect. 4.5, and for their comments on a first version of this manuscript; R. Capuzzo-Dolcetta, B. Famaey, V. Hill, R. Ibata, D. Katz, N. Martin, A. Mastrobuono-Battisti, M. D. Lehnert, Y. Revaz, and O. N. Snaith, for stimulating discussions and remarks. This work was granted access to the HPC resources of TGCC and CINES under the allocations 2014-040507 and 2016-040507 made by GENCI. Finally, we wish to thank the referees, for their prompt reports and for their comments which helped improving the clarity of this paper.

References

Appendix A: The “isolated” Milky Way

To understand how one or several mergers can impact the distribution of stars in integrals-of-motion and kinematic spaces, it is interesting to compare with the distribution one would obtain for a Milky Way-type galaxy evolved in isolation. To achieve this aim, we ran a simulation of the Milky Way-type galaxy, evolved in isolation for 5 Gyr and which initially has the same internal parameters as those adopted for the minor merger simulations (see Table 1). In Figs. A.1–A.3 the corresponding distribution of stars in the E − L_z space, in the L_⊥ − L_z space, and in the R_apo − R_peri space are shown. Comparing these figures with those presented in Sect. 3, it is clear that interactions tend to redistribute stars over a much larger portion of all the spaces considered. In the E − L_z plane, halo stars in the isolated galaxy, which constitute the tail of the thick disc distribution, are few and tend to redistribute over a thin region in this space. In the L_⊥ − L_z space, halo stars in solar vicinity volumes redistribute over a considerably less extended region of the space than that occupied in the case of one or several accretion events (cfr Figs. A.2 and 7). A similar conclusion can be reached for the R_apo − R_peri space: looking at the distribution of all stellar particles in the simulation, three main structures are observed, corresponding to thin disc, intermediate disc and thick disc stars (see left panel, Fig. A.3). When the analysis is restricted to the solar vicinity volume (middle and right panels, Fig. A.3), one sees how much thinner the distribution is with respect to those shown in Figs. C.6 and C.7. Note in particular the presence of stars belonging to the stellar bar, which show up as an over-density at high eccentricity in this plane (R_peri< 2 kpc, R_apo< 10 kpc).

Appendix B: The spatial distribution of clumps in the E – L_z space

In this section we show the spatial distribution of some of the stellar clumps observed in the E − L_z space. To this end, in the top panels of Figs. B.1–B.3 we report the distribution in this space for all stars, in-situ stars only and accreted stars only, respectively, in the case of the 2 × (1:10) merger. From this distribution we have selected by eye, four clumps, which correspond to local maxima in the E − L_z space and indicated by the symbols #1, #2, #3 and #4 in all the figures. For all stars belonging to each of these clumps, we report in the second, third, fourth and fifth rows of Figs. B.1–B.3 (in order of decreasing energy) their spatial distribution in the x − y and in the x − z planes.

As expected, for all populations (in-situ, accreted and the two combined), the spatial distribution of stars in a given clump depends on its energy, with stars in the less densely populated clumps having a more extended radial and vertical extent. Also, stars belonging to less densely populated clumps show a less homogeneous (and unrelaxed) spatial distribution than stars belonging to clumps of lower energy, no m atter what their origin, in-situ or accreted. Finally we emphasize that both in-situ and accreted stars show a very disturbed spatial distribution in the outer disc: stellar streams, plumes, are indeed visible in both populations, as is evident by comparing the x − z distribution in Figs. B.2 and B.3.

Appendix C: Additional figures

All Tables

All Figures