<figure><img src="image_1.png" /></figure>

# In Disguise or Out of Reach: First Clues about In Situ and Accreted Stars in the Stellar Halo of the Milky Way from Gaia DR2

M. Haywood<fnref n="1" />, P. Di Matteo<fnref n="1" />, M. D. Lehnert<fnref n="2" />, O. Snaith<fnref n="3" />, S. Khoperskov<fnref n="1" />, and A. Gómez<fnref n="1" />

Received 2018 May 8; revised 2018 July 5; accepted 2018 July 5; published 2018 August 16

## Abstract

We investigate the nature of the double color–magnitude sequence observed in the Gaia DR2 HR diagram of stars with high transverse velocities. The stars in the reddest-color sequence are likely dominated by the dynamically hot tail of the thick disk population. Information from Nissen & Schuster and from the APOGEE survey suggests that stars in the blue-color sequence have elemental abundance patterns that can be explained by this population having a relatively low star formation efficiency during its formation. In dynamical and orbital spaces, such as the “Toomre diagram,” the two sequences show a significant overlap, but with a tendency for stars on the blue-color sequence to dominate regions with no or retrograde rotation and high total orbital energy. In the plane defined by the maximal vertical excursion of the orbits versus their apocenters, stars of both sequences redistribute into discrete wedges. We conclude that stars that are typically assigned to the halo in the solar vicinity are actually both accreted stars lying along the blue sequence in the HR diagram, and the low rotational velocity tail of the old Galactic disk, possibly dynamically heated by past accretion events. Our results imply that a halo population formed in situ and responsible for the early chemical enrichment prior to the formation of the thick disk has yet to be robustly identified, and that what has been defined as the stars of the in situ stellar halo of the Galaxy may in fact be fossil records of its last significant merger.

**Key words:** Galaxy: evolution – Galaxy: halo – Galaxy: kinematics and dynamics

## 1. Introduction

Over the past three decades, a consensus has developed that the stellar halo of the Milky Way is composed of two populations of stars—those that were born in other galaxies and accreted and those that were born in situ during the early evolution of the Milky Way (e.g., Searle & Zinn 1978; Sommer-Larsen & Zhen 1990; Carollo et al. 2007). Although these two components have remained challenging to characterize, we now have a widely accepted picture whereby the in situ halo population is older, is more metal-rich, dominates the stellar density within ~15 kpc of the Galactic center, and has a slightly enhanced mean rotation rate compared to the accreted halo population (e.g., Carollo et al. 2007). While a considerable amount of effort is currently being expended to investigate stellar streams that are expected to fill the outer stellar halo, the inner halo population is still lacking a proper characterization, and its evolutionary connection with the thick disk is essentially unknown.

The second data release (Gaia Collaboration et al. 2018b) of the European Space Agency’s Gaia mission (Gaia Collaboration et al. 2016) provides superb astrometric parameters, radial velocities, and photometry for a large number of stars. First results have already detected a number of streams and kinematic groups (Koppelman et al. 2018; Malhan et al. 2018). Inspection of the Gaia HR diagram (hereafter HRD) of stars with high total or tangential velocities has shown two parallel color sequences (Gaia Collaboration et al. 2018a), which were attributed to stars in the thick disk and stellar halo. Here, we concentrate on understanding the origin of these two sequences and in particular the nature of the blue-color sequence. The next section describes our selection from the Gaia archive, Section 2 analyzes the sample of stars from Nissen &

Schuster (2010) specifically with respect to where stars in this sample lie in these two sequences and thus bringing insights into the characteristics of the blue-color sequence. In Section 2.3, we present their main kinematic and orbital properties and we discuss our results and summarize our main conclusions in Section 4.

## 2. Clues from Nissen & Schuster (2010)

### 2.1. Data

We select stars in the Gaia archive having a tangential velocity greater than 200 km s⁻<fnref n="1" />, parallax errors σ<sub>π</sub>/π < 0.1, π > 1 mas, G < 17, and the filters described in Gaia Collaboration et al. (2018a). This brings us 77,107 stars displayed in the HRD of Figure 1 (left panel). This figure shows that the HRD separates in two sequences, as in Gaia Collaboration et al. (2018a), and which we designate below as the blue and the red sequence (BS and RS, respectively).

In order to clean the HRD, we select stars with various levels of tolerance on the interstellar reddening. We adopt the interstellar reddening estimates from the map of Lallement et al. (2018), which is well adapted for stars nearer than 1 kpc. Figure 1 (middle and right panels) shows the HRD obtained by selecting stars that have levels of reddening below 0.025 and 0.015 mag, providing, respectively, 28210 and 12620 objects.

### 2.2. Chemical Abundances and HR Diagram

The sample of Nissen & Schuster (2010; hereafter NS) in three different metallicity intervals is overplotted to the Gaia HRD in Figure 2, together with isochrones from the PARSEC library (Marigo et al. 2017) at Z = (0.006, 0.0024, 0.0048) and ages (11, 11.5, 12.5) Gyr. Note that the

<footnote n="1">¹GEPI, Observatoire de Paris, PSL Research University, CNRS, Sorbonne Paris Cité, 5 place Jules Janssen, F-92190 Meudon, France; misha.haywood@obspm.fr</footnote>

<footnote n="2">²Sorbonne Université, CNRS UMR7095, Institut d'Astrophysique de Paris, 98 bis bd Arago, F-75014 Paris, France</footnote>

<footnote n="3">³School of Physics, Korea Institute for Advanced Study, 85 Hoegiro, Dongdaemun-gu, Seoul 02455, Republic Of Korea</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** HRD of stars selected in Section 2.1, with different tolerances in the reddening. The two, blue- and red-color sequences, are well separated.</figcaption></figure>

Gaia magnitudes were derived using the V, I magnitudes and color transformations provided by Evans et al. (2018).

This figure shows that the red and blue sequences are dominated by metal-rich ([Fe/H] > −0.8) and metal-poor ([Fe/H] < −1.1) stars (as already anticipated in Gaia Collaboration et al. 2018a). Note that the spread in metallicity in each sequence is significant, ~0.6 dex. A cross-match between APOGEE (Majewski et al. 2017) and our Gaia sample yields 226 stars, showing that there is a marked dip in the metallicity distribution of stars at [Fe/H] ~ −1.0 (Gaia Collaboration et al. 2018a). This dip is also noticeable in the data of NS (their Figure 1). The middle plot of Figure 2 shows that stars with −1.1 < [Fe/H] < −0.8, which brackets the dip in metallicity, fall between the two sequences of the HRD. This implies that the clear separation between the two sequences is the conspicuous consequence of the dip in the metallicity distribution function (MDF). Moreover, the middle plot shows that, at similar metallicities, low-α stars tend to be bluer than high-α ones, suggesting that they must be slightly younger, supporting similar results from Schuster et al. (2012). On the contrary, the bottom panel of Figure 2 shows that low- and high-α stars with [Fe/H] < −1.1 from the NS sample can both be found along the BS stars. This is surprising, because—given the tightness of the blue-color sequence—if the low- and high-α stars at [Fe/H] < −1.1 belonged to two different chemical evolution sequences—as they do at −1.1 < [Fe/H] < −0.8— we would expect a difference in age that would reflect in the HRD, as is seen for their more metal-rich counterparts (Figure 2, middle plot). Is it thus possible that all the stars with [Fe/H] < −1.1—both high- or low-α stars, because the distinction becomes less obvious below this metallicity—are the same population of stars, and are not causally related to the thick disk? Could they instead be causally related to the low-α sequence at higher metallicities, which is clearly distinct from the thick disk?

By analyzing the [Fe/H]–[Mg/Fe] distribution of APOGEE stars, Hayes et al. (2018) have shown that the separation between the high- and low-α stars is not horizontal—it is not a separation in α-abundances—but that high-α stars at [Fe/H] ≲ −1.1 and low-α stars at higher metallicities form a unique chemical sequence. The findings of Hayes et al. (2018) strongly support our suggestion. Figure 3⁴ shows the [Fe/H]–

[α/Fe] distribution of the sample of Nissen & Schuster (2010), with the high-α, low-metallicity stars shown as a blue square. Indeed, these objects are all within the BS. It is therefore reasonable to suggest that the low-metallicity ([Fe/H] < −1.1) stars of NS and the low-α sample at higher metallicity are forming a unique abundance sequence. The gap in the MDF at [Fe/H] ~ −1 could simply be a reflection of the transition between two populations of stars: a population whose origin needs to be determined (see below) at [Fe/H] < −1, and the thick disk above this limit. Hayes et al. (2018) also mention that the transition between the two sequences in the [Fe/H]–[α/Fe] plane is more pronounced at this metallicity (see also Bonaca et al. 2017).

## 2.3. Kinematics and Orbits

Of the sample of 28,210 stars with $E(B-V) < 0.025$ mag, 1973 stars have full 3D velocity information in Gaia DR2. In the following, we have assumed an in-plane distance of the Sun from the Galactic center $R_⊙ = 8.34$ kpc following Reid et al. (2014), a height of the Sun above the Galactic plane $z_⊙ = 27$ pc (Chen et al. 2001), a velocity for the Local Standard of Rest $V_{LSR} = 240$ km s⁻¹ (Reid et al. 2014), and a peculiar velocity of the Sun with respect to the LSR, $U_⊙ = 11.1$ km s⁻¹, $V_⊙ = 12.24$ km s⁻¹, $V_⊙ = 7.25$ km s⁻¹, following Schönrich et al. (2010).

Figure 4 shows the HR diagram of the subsample of stars with 3D kinematics, and how we separate them using the isochrone.

The distributions of our stars in kinematic spaces is given in Figure 5. It shows two overlapping structures: a vertical plume that extends to high energies and was already noted in Gaia Collaboration et al. (2018c, their Figure 5), and is composed of stars on the blue sequence, as noted by Koppelman et al. (2018), and the high energy tail of the disk, which extends to higher rotation, and contains stars on the red sequence. The two sequences overlap in the Toomre diagram; however, the majority of stars on the RS, ≈75%, shows prograde rotation and lies in the region of the diagram that is compatible with the slow rotating tail of the thick disk (Figure 5). The remaining fraction, 1/4 of red-color sequence stars, however, has a retrograde orbit, lagging the LSR by ~400 km s⁻¹ or less. Except for 13 stars with very significant retrograde motions, $V < -500$ km s⁻¹, that we will discuss in the following, stars on the BS are distributed over the same area as the RS stars in the Toomre diagram, but in different proportions. In particular, they dominate the extended vertical plume at null or retrograde rotation, with about 65% of BS stars lying within the region with $-350$ km s⁻¹ $≤ V ≤ -200$ km s⁻¹, while 26% have $V > -200$ km s⁻¹ and the remaining fraction $V < -350$ km s⁻¹.

<footnote n="4">⁴ In the first version of this article, G 81-02 was found to fall blueward to the blue sequence, due to a probably erroneous extinction estimate. P. Nissen informed us that his own estimate from uvby photometry indicates a reddening of $E(G_{BP}-G_{RP}) = 0.043$ mag, 5 times less than ours. This seems to be much more consistent with the metallicity of this star (−0.69 dex), and has been used here.</footnote>
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 2.** Overplotted on the Gaia HRD, stars from Nissen & Schuster (2010) with [Fe/H] > −0.8 (top panel), −1.1 < [Fe/H] < −0.8 (middle panel), [Fe/H] < −1.1 (bottom panel). Thick disk and high-α halo stars from NS are represented with red circles, low-α halo stars from NS with blue triangles. The isochrones from the PARSEC library have (top to bottom) (0.006 dex, 11 Gyr), (0.0048 dex, 11.5 Gyr), and (0.0024 dex, 12 Gyr).</figcaption></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 3.** Top: the low- and high-α stars from Nissen & Schuster (2010, triangles and squares). The high-α stars at [Fe/H] < −1.2 plotted as blue squares may be assumed to be part of the same chemical sequence (Hayes et al. 2018). Bottom: the resulting classification is more consistent with that belonging to the two sequences in the HRD, see Figure 2.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 4.** Gaia HR diagram of the stars with 3D kinematics in our sample. We separate our sample into two groups (red and blue sequence) according to the position of each star relative to the isochrone. The isochrone is the same as the one in Figure 2, middle plot.</figcaption></figure>

Remarkably, stars of the BS with prograde orbits clearly overlap with those of RS stars.

In order to understand how these results are sensitive to the adopted classification on the red and blue sequences, we studied how the Toomre diagram changes when adopting a stricter separation. Stars have been selected, respectively,
<figure><img src="image_9.png" /><figcaption>**Figure 5.** From top to bottom: Toomre diagram, $E-L_z$, $L_z-L_{perp}$ and $z_{max}-R_{2D,max}$ planes for all stars in our sample with full 3D velocities in Gaia DR2. The dashed line in the top plot indicates $V = -V_{LSR}$. In the bottom panel, only stars with $R_{2D,max} \lesssim 40$ kpc have been plotted. In all panels, stars on the blue sequence are shown as filled blue dots, stars on the red sequence are shown as empty red circles.</figcaption></figure>

+0.05 redder and −0.05 bluer than the isochrone. The result is shown in Figure 6 and shows that the effect is slightly different for the two populations: the contours of the distribution of stars on the blue sequence remain basically the same, with even the stretch of stars at strongly retrograde orbits being still present in the selection. On the contrary, several stars of the red sequence on high energy orbits have been removed, suggesting that the vertical plume, at least at the highest energies, could be almost exclusively populated by stars on the blue sequence. While a stricter separation preferentially removes stars on high energy orbits from the red sequence, we note that apart from this effect, the main characteristics of the two distributions are essentially unchanged using the two different selections.

The same main structures visible in the Toomre diagram are also found in the $E-L_z$ plane (Figure 5), where the energy $E$ is the total energy of a star, defined as the sum of its kinetic and gravitational potential energies. We assume an Allen & Santillan (1991) (hereafter AS) Galactic mass model to estimate the latter. In this plane, we note that the RS stars around $L_z = 0$ are predominately on very bound orbits (i.e., low $E$). In the $L_z-L_{perp}$ plane, where $L_{perp} = \sqrt{L_x^2 + L_y^2}$, there is significant overlap between the BS and RS within this space, except in two regions: (1) the region where the most extreme counter-rotating stars lie, which is composed of stars exclusively from the blue-color sequence; and (2) an over-dense region roughly centered at $(L_z, L_{perp}) \approx (1500, 2200)$ kpc km s⁻¹, which is very close to the region occupied by the Helmi stream ($(L_z, L_{perp}) \approx (1000, 2000)$ kpc km s⁻¹, see Helmi et al. 1999), which is also exclusively made of stars on the BS.

Finally, for all stars with full 3D velocity information, we have reconstructed their orbital parameters by integrating their orbits over the last 6 Gyr, using four different Galaxy potentials: the axisymmetric potential of Allen & Santillan (1991), the two axisymmetric potentials, including a thick disk component, in Pouliasis et al. (2017), and the MWpotential2014 from Bovy (2015). In Figure 5, we show their distribution in the $z_{max}-R_{2D,max}$ plane, where $z_{max}$ is the maximum height that stars reach above or below the Galactic plane and, $R_{2D,max}$, is the apocenter of their orbit projected on the Galactic plane. For clarity, we show in this plot only the orbital parameters that have been derived by integrating the orbits in the Allen & Santillan (1991) potential, and below we discuss the robustness of this analysis when the other three potentials are used. Predominately, the stars on the RS have $R_{2D,max}$ within 20 kpc of the Galactic center, while the orbits of stars on the BS can extend much further out than 20 kpc. A striking feature is that the stars in our sample are not homogeneously distributed in this plane, but define three distinct diagonal “wedges,” with $z_{max}$ increasing with $R_{2D,max}$.

Remarkably, one of these three patterns defines stars confined in a relatively thin and flattened distribution, with a lack of stars between $2 \le z_{max} \le 4$ kpc. This lack of stars with $z_{max}$ in this range makes this region distinct from the rest of the sample. Of the 13 stars with very significant retrograde motions, $V < -500$ km s⁻¹, 10 are found in this thin flattened disk. Note that the presence of two groups of halo stars, a first with a flattened distribution and a second with more vertically extended orbit, was already noted by Schuster et al. (2012, see their Figure 8).
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Top: Toomre diagram shown separately for stars on the blue and red sequences, with the same separation as in Figure 5. Bottom: the separation is now done by selecting stars bluer by -0.05 than the isochrone and redder by +0.05 mag for the stars in red, so the two samples are separated by 0.1 mag in color, decreasing possible contamination between the two. See the text for details.</figcaption></figure>

In order to assess if the structures in the (RD2max, Zmax) plane are robust to (1) the uncertainties in the parallax, proper motions, and radial velocity; and (2) the uncertainties in the mass distribution in the Galaxy, we generated 10 random sets of these parameters for each star around the observed values, using the $1\sigma$ uncertainties given in the Gaia DR2 catalog. The orbits were then integrated from the *U*, *V*, *W* velocities in four different Galactic potentials described above. From the orbital parameters generated in all four potentials, we derived the angles defined by the arctan(Zmax/R2Dmax), and plot the resulting distributions in Figure 7.

The potential of Model I in Pouliasis et al. (2017) is similar to the Allen & Santillan (1991) potential, but includes a massive thick disk; however, since the overall mass distribution is fitted to the same rotation curve, the overall potential is not very different, hence the distributions of Figure 7 (top) are very similar. On the contrary, Model II in Pouliasis et al. (2017) was fitted to a different rotation curve (see Pouliasis et al. 2017, for details), having a more massive dark matter halo, and thus its distribution differs significantly from the previous two potentials, but is similar to the distribution obtained with the MWpotential2014 from Bovy (2015).

A common feature seen in all four distributions is the dip between 10° and 25° (Figure 7). It is wider in the first two potentials, but is significant in all resulting distributions. The gap between the intermediate and highest wedge is visible in

only one potential, but there seems to be nothing significant in the other three potentials. We conclude that (only) the dip at 10°–25° is significant.

## 2.4. The Effect of the Kinematic Selection

In this study, we selected stars having tangential velocities higher than 200 km s⁻¹. Does the cut on tangential velocities affect our conclusions? One possible way to answer this question is to compare our results with samples that were not kinematically defined. To do so, we look at the data set of Chiba & Beers (2000), which has been designed to avoid kinematic bias. Figure 8 (top) shows the (R2Dmax, Zmax) plane adopting the orbital parameters published in Chiba & Beers (2000), which shows that the structures were already apparent in the (Zmax, R2Dmax) plane. Figure 8 (bottom) shows the same distribution but with orbital parameters recalculated in the AS potential using the Gaia DR2 astrometric parameters. The separation between the lower and intermediate wedges is clear. The stars from Schuster et al. (2012), comprising halo stars with total velocities greater than 180 km s⁻¹ but also thick disk stars (see Figure 3 in Nissen & Schuster 2010), are also overplotted (black triangles), with most objects populating the low and intermediate wedges, with a clear separation between the two. Note that the orbital parameters of the Schuster et al. (2012) stars shown in Figure 8
<figure><img src="image_11.png" /><figcaption>**Figure 7.** Distribution of angles (arctan(Zmax/R2Dmax)) obtained by calculating 10 orbits for each star in four different Galactic potentials. The dip in the distribution (between 10° and 25°) between the lower and intermediate wedges is found in all four potentials. The dip between the intermediate and upper wedges is not.</figcaption></figure>

(bottom) correspond to those derived in an asymmetric barred potential (see Schuster et al. 2012 for details).

### 3. The APOGEE-Gaia Sample

In order to further explore the picture presented above, there are several limitations of our analysis so far that we can overcome. Since our sample is limited to distances where stars have $\pi > 1$ mas, our sample contains few giants, which are the objects with the best S/N in APOGEE. Therefore, in order to include more giants, we selected in the Gaia DR2 catalog stars having $\pi > 0.3$ mas, $\sigma_{\pi}/\pi > 0.1$ and $G < 15$. The corresponding color-magnitude diagram is shown in Figure 9, not corrected for extinction and reddening (top) and corrected (bottom). Both sequences are well separated and densely populated and contain considerably more giants, as expected. Cross-matching this Gaia sample with APOGEE yields 950 stars with abundance estimates. Figure 10 gives the (R2Dmax, Zmax), (E,Lz) distributions and Toomre diagram of these 950 stars, as calculated in the AS potential. As expected, the new sample closely follows the distribution of Figure 5.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Top: the (R2Dmax, Zmax) distribution of stars from Chiba & Beers (2000) with the orbital parameters determined by the authors. The same sample with the orbital parameters recalculated with the AS potential and Gaia astrometric parameters and radial velocity, together with the sample of Schuster et al. (2012; black triangles), with their parameters derived from orbits calculated in the barred potential.</figcaption></figure>

Based on what we learned in the previous section, we now select objects with R2Dmax > 20 kpc on the three “wedges” colored in blue, red, and green, assuming these are the most likely to have been accreted, with the aim to see how they distribute in the [Fe/H]–[Mg/Fe] plane. In the lower altitude wedge, stars are represented by two symbols according to the value of their pericenter: below or above 2 kpc in blue or cyan respectively.

The [Fe/H]–[Mg/Fe] distribution of these different sub-samples is given in Figure 11. It shows that stars with large apocenters (high energy orbits) clearly result in forming a low star formation efficiency sequence, which is typical of a relatively massive dwarf galaxy, extending from [Fe/H] ~ -0.7 to [Fe/H] ~ -2.0.

Stars in the low altitude wedge are a mix of different origins: those colored in blue also belong to the low star formation efficiency sequence, and they differ from the others (colored in cyan) by having pericenters smaller than 2 kpc and which bring them near to the central regions. The objects colored in cyan have larger pericenters, typically between 3 and 7 kpc, and are mostly thick disk stars (Figure 11). They represent the tail of this population scattered at large galactocentric distances.

Figure 10 shows as magenta dots stars that are chemically defined as thick disk objects in Figure 11. It confirms that bona
<figure><img src="image_13.png" /><figcaption>**Figure 9.** *Gaia* color–magnitude diagram for the second sample described in Section 3, not corrected for extinction (top) and corrected (bottom). The black dots on the bottom plot are the 950 stars found in the cross-match between the *Gaia* sample and APOGEE.</figcaption></figure>

fide thick disk stars have a rotation that decreases to zero, and some are counter-rotating, up to relatively high metallicities ($[Fe/H] > -0.5$).

Finally, Figure 12 shows the stars with large apocenters and those belonging to the thick disk overplotted on the color–magnitude diagrams of the *Gaia* sample. All stars have been corrected for extinction. The figure confirms the previous analysis by showing that the blue sequence is dominated by stars that are on the accreted sequence, while the red sequence is dominated by thick disk objects, including those that have large apocenters but clearly belong to the chemically defined thick disk population.

<figure><img src="image_14.png" /><figcaption>**Figure 10.** Top: (R2Dmax, Zmax) distribution of stars in the APOGEE-*Gaia* cross-match. Stars belonging to the three wedges with R2Dmax > 20 kpc are colored in green, red, and blue. The stars of the low altitude wedge are plotted in two different colors according to the pericenter: blue for stars with R2Dmin < 2 kpc and cyan for stars with R2Dmin > 2 kpc. In magenta are the thick disk stars selected to have $[Fe/H] > -1$ and $[Mg/Fe] > 0.2$ in the next figure. Middle and bottom plots show the (E,Lz) distribution and Toomre diagram for these stars.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Figure 11.** $[Fe/H]–[Mg/Fe]$ distribution of stars of the APOGEE-*Gaia* sample. The stars selected on the basis of their R2Dmax (blue, red, and green) are the same as in Figure 10. Thick disk stars are represented by magenta points.</figcaption></figure>

## 4. Discussion and Conclusions

Our results come in the context of a long history of halo population studies. The dichotomy observed in the *Gaia* HRD
<figure><img src="image_16.png" /><figcaption>**Figure 12.** *Gaia* color–magnitude diagram for the APOGEE stars selected as described in Section 3. Symbols are the same as those in Figure 10. Magnitudes have been corrected for extinction.</figcaption></figure>

has been known for several decades as due to the halo and thick disk. Gilmore et al. (1985) found that star counts at the pole are dominated by two populations with main-sequence turn-offs dominated by stars at $[Fe/H]$-1.5, and the other at $[Fe/H] = -0.7$. This is also what was determined in the analysis of Sloan data (Ivezić et al. 2008). More recently, Jofré & Weiss (2011) studied the main-sequence turn-off of old populations using SDSS data and determined the age of the inner Galactic halo. The flattened distribution found here, which reflects a strongly anisotropic velocity distribution, is reminiscent of the findings of Sommer-Larsen & Zhen (1990) and echoes the analysis of Chiba & Beers (2000). After their results, the Galactic halo sampled at the solar vicinity has been seen as an inner, flattened, and in situ counterpart of a more extended, spherical, and accreted halo, a view that has been supported by other, more recent studies (see Carollo et al. 2007, 2010; Beers et al. 2012), but it must be noted that Brook et al. (2003) already proposed that eccentric stars in the sample of Chiba & Beers (2000) could have been the result of an accretion. The main result of our study is that the inner halo is probably mainly composed of accreted material left by the last significant merger in our Galaxy.

Our results confirm that the red sequence must be dominated by thick disk stars with metallicities between −0.4 and −1, as found by Gaia Collaboration et al. (2018a). The majority of RS stars are on prograde orbits. However, an important, new result is that a non-negligible fraction of RS stars are on retrograde orbits. This can be at least partially attributed, in Figure 5, to an imperfect separation of stars of the BS and RS in the HRD—that is, part of the RS stars with very low or retrograde tangential velocities may in fact be BS stars that have been incorrectly assigned to the RS. Importantly, the analysis of the APOGEE sample clearly shows that some stars defined chemically as thick disk objects have retrograde orbits. It is known that in situ stellar disks, dynamically heated by one or several satellite accretion events, will also have counter-

rotating stars (see Qu et al. 2011; Jean-Baptiste et al. 2017). Hence, finding some nonrotating or counter-rotating stars among an in situ old disk population is not surprising. Although they have chemical characteristics typical of the thick disk, RS stars have the kinematics of halo stars (see also Bonaca et al. 2017; Posti et al. 2018). The most likely origin of the red-color sequence is thus related to the old Galactic disk, partially heated during some past accretion event(s).

Gaia Collaboration et al. (2018a) suggest that the blue-color sequence may be related to the in situ halo. We find here that this sequence is most probably dominated by accreted objects, as also suggested by Koppelman et al. (2018).

Low-metallicity, $[Fe/H] < -1.1$, high-$\alpha$ stars from Nissen & Schuster (2010) and from the APOGEE sample studied here are indeed on the BS. Our analysis shows that low-metallicity, low-$\alpha$ stars are also on this sequence (Section 2). Based on the fact that they belong to the same narrow sequence in the HRD and on the results of Hayes et al. (2018), we conjecture that the high-$\alpha$ stars with $[Fe/H] < -1.1$ could be causally connected to the higher metallicity, low-$\alpha$ stars forming a unique chemical evolution sequence. This is confirmed by the chemical characteristics of the stars with the highest orbital energies in our APOGEE sample: these objects form a sequences that extends from about $[Fe/H] = -0.5$ to at least $-2.0$ and from low- to high-$\alpha$ abundance. Most of these objects are on the blue sequence in the *Gaia* HRD.

If these stars are on the same evolutionary sequence, they must have been formed in an environment characterized by a low star formation efficiency (i.e., long star formation timescale), much lower than that of the thick disk. For this sequence, Supernovae-Ia must have started enriching the interstellar medium in Fe at lower metallicities compared to what is observed for thick disk stars of the Milky Way. The low-$\alpha$ stars in the sample of Nissen & Schuster (2010) would then only represent the “tip of the iceberg” of the much larger population blue-color sequence stars studied by us and APOGEE (Hayes et al. 2018). The most natural origin of the blue-color sequence is that it is dominated by stars accreted from a satellite.

The stars that most likely were accreted dominate the distribution of (Zmax,Rmax) at the largest apocenters, but it is difficult to determine the fraction they represent at lower R2Dmax. Clearly, the low-$\alpha$ sequence is not only populated by objects with large apocenters. Interestingly, the sample at metallicities lower than −2 have more limited excursions in R and Z, with typical orbits being limited to within 15 and 10 kpc. This peculiarity has also been noted by Schuster et al. (2012).

What might be the origin of these accreted stars? Belokurov et al. (2018) find a strong orbital anisotropy for stars with metallicities above −1.7 and low anisotropy for stars below this limit. They suggest that the majority of the halo stars within 30 kpc are remnants of a massive satellite accreted during the formation of the Galactic disk between about 8 and 11 Gyr ago. The analysis of Gaia DR2 data may support this scenario for several reasons: (1) BS stars seem to constitute a significant fraction of all stars with high transverse velocities; (2) we show that they could be on a chemical evolutionary track that is less $\alpha$-enriched than disk stars at the same metallicity and thus are compatible with an accreted population; (3) the kinematic properties of this BS—in particular, the high fraction of stars in retrograde orbits, a fraction of which are confined within a flattened and extended disk—and more generally the discrete wedges in the $z_{max}-R_{2D,max}$ plane—are all reminiscent of some impulsive heating of the early Galactic
disk related to some accretion event(s). In this respect, the gap in the distribution of $z_{max}$ values may mark the transition from an early phase of significant stellar accretion in the Galaxy to a more quiescent phase.

Can all stars in the blue-color sequence be attributed to a unique satellite accretion event, or is it possible that several satellites contributed to make it? From the tightness of the BS in the HRD, it is difficult to conceive that this sequence consists of populations formed in several satellite galaxies, unless their chemical and age properties at the time of their accretion were remarkably similar. Also the kinematics of stars in this sequence may be compatible with a unique—and relatively massive—merger. Figure C.3 in Jean-Baptiste et al. (2017), for example, shows that among several satellites accreted onto a Milky Way-type galaxy, stars originating in one of them (satellite #3 in that plot) are distributed in the $E-L_z$ space in a way qualitatively similar to the BS stars we have observed, with stars having both prograde and retrograde orbits, and a significant plume of stars with $L_z$ centered around $L_z = 0$, and extending vertically to high energies. These simulations support the notion that the blue-color sequence is (at least partially) the remnant of a significant accretion event in the early history of the Milky Way. Note that Nissen & Schuster (2010) and Schuster et al. (2012) have evoked the possibility that the low-$\alpha$ stars in their work may originate from $\omega$ Cen, and at this stage we cannot reject or confirm this suggestion for the origin of the blue sequence. Koppelman et al. (2018) also note that previous studies have associated this region to possible debris from $\omega$ Cen. However, we note that $\omega$ Cen seems to have a peculiar barium abundance (Majewski et al. 2012) that does not appear to be compatible with that of the low-$\alpha$ stars discussed here (see Nissen & Schuster 2011).

While our work represents only a first exploration of the stellar halo in *Gaia* DR2, it is clear that this mission, together with large spectroscopic surveys, is reshaping the boundaries that we had for decades assigned to the various stellar populations of the Milky Way. Our results suggest that what has been defined as the stars of the in situ stellar halo of the Galaxy may be in fact fossil records of its last significant merger. Stars kinematically defined to belong to the inner halo comprise this possible accreted population and a sizable fraction of more metal-rich, $[Fe/H] \gtrsim -1$, stars that are possibly the vestiges of the early disk of the Galaxy after it was heated by one or more merger event. But where is the parent population of the thick disk? Is it possible that for this progenitor, the in situ, chemically defined halo is lurking in the blue sequence, but is under-represented in the volume probed by our study? So the question remains, is the in situ halo stellar population expected in galaxy evolution models disguised among the blue-color sequence or is it still beyond our reach?

We would like to thank Poul Nissen and Timothy Beers for valuable comments, the referee for several suggestions that improved this work, and David Katz and Nicolas Leclerc for helping with the use of the *Gaia* Archive. The Agence Nationale de la Recherche (ANR) is acknowledged for its financial support through the MOD4Gaia project (ANR-15-CE31-0007, P. I.: P. Di Matteo), also for providing the postdoctoral grant for Sergey Khoperskov. This work has made use of data from the European Space Agency (ESA) mission *Gaia* (https://www.cosmos.esa.int/gaia), processed by the *Gaia* Data Processing and Analysis Consortium (DPAC, https://www.cosmos.esa.int/web/gaia/dpac/consortium).

Funding for DPAC has been provided by national institutions, in particular, the institutions participating in the *Gaia* Multilateral Agreement. This research has made use of the SIMBAD database, operated at CDS, Strasbourg, France. Funding for the Sloan Digital Sky Survey IV has been provided by the Alfred P. Sloan Foundation, the U.S. Department of Energy Office of Science, and the Participating Institutions. SDSS-IV acknowledges support and resources from the Center for High-Performance Computing at the University of Utah. The SDSS website is www.sdss.org. SDSS-IV is managed by the Astrophysical Research Consortium for the Participating Institutions of the SDSS Collaboration, including the Brazilian Participation Group, the Carnegie Institution for Science, Carnegie Mellon University, the Chilean Participation Group, the French Participation Group, Harvard-Smithsonian Center for Astrophysics, Instituto de Astrofísica de Canarias, The Johns Hopkins University, Kavli Institute for the Physics and Mathematics of the Universe (IPMU)/University of Tokyo, Lawrence Berkeley National Laboratory, Leibniz Institut für Astrophysik Potsdam (AIP), Max-Planck-Institut für Astro-nomie (MPIA Heidelberg), Max-Planck-Institut für Astrophysik (MPA Garching), Max-Planck-Institut für Extraterrestrische Physik (MPE), National Astronomical Observatories of China, New Mexico State University, New York University, University of Notre Dame, Observatório Nacional/MCTI, The Ohio State University, Pennsylvania State University, Shanghai Astronom-ical Observatory, United Kingdom Participation Group, Uni-versidad Nacional Autónoma de México, University of Arizona, University of Colorado Boulder, University of Oxford, Uni-versity of Portsmouth, University of Utah, University of Virginia, University of Washington, University of Wisconsin, Vanderbilt University, and Yale University.

*Note.* A month after this article in its first version was submitted to this journal and posted on arXiv, Helmi et al. posted an article that confirms some of our results.

## ORCID iDs

M. Haywood https://orcid.org/0000-0003-0434-0400

P. Di Matteo https://orcid.org/0000-0002-5213-4807

M. D. Lehnert https://orcid.org/0000-0003-1939-5885

## References

* Allen, C., & Santillan, A. 1991, RMxAA, 22, 255

* Beers, T. C., Carollo, D., Ivezić, Ž., et al. 2012, ApJ, 746, 34

* Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, MNRAS, 478, 611

* Bonaca, A., Conroy, C., Wetzel, A., Hopkins, P. F., & Kereš, D. 2017, ApJ, 845, 101

* Bovy, J. 2015, ApJS, 216, 29

* Brook, C. B., Kawata, D., Gibson, B. K., & Flynn, C. 2003, ApJL, 585, L125

* Carollo, D., Beers, T. C., Chiba, M., et al. 2010, ApJ, 712, 692

* Carollo, D., Beers, T. C., Lee, Y. S., et al. 2007, Natur, 450, 1020

* Chen, B., Stoughton, C., Smith, J. A., et al. 2001, ApJ, 553, 184

* Chiba, M., & Beers, T. C. 2000, AJ, 119, 2843

* Evans, D. W., Riello, M., De Angeli, F., et al. 2018, arXiv:1804.09368

* Gaia Collaboration, Babusiaux, C., van Leeuwen, F., et al. 2018a, arXiv:1804.09378

* Gaia Collaboration, Brown, A. G. A., Vallenari, A., et al. 2018b, arXiv:1804.09365

* Gaia Collaboration, Katz, D., Antoja, T., et al. 2018c, arXiv:1804.09380

* Gaia Collaboration, Prusti, T., de Bruijne, J. H. J., et al. 2016, A&A, 595, A1

* Gilmore, G., Reid, N., & Hewett, P. 1985, MNRAS, 213, 257

* Hayes, C. R., Majewski, S. R., Shetrone, M., et al. 2018, ApJ, 852, 49
Haywood, M., Di Matteo, P., Lehnert, M., Katz, D., & Gómez, A. 2013, A&A, 560, A109

Haywood, M., Di Matteo, P., Snaith, O., & Calamida, A. 2016, A&A, 593, A82

Helmi, A., White, S. D. M., de Zeeuw, P. T., & Zhao, H. 1999, Natur, 402, 53

Ivezić, Ž., Sesar, B., Jurić, M., et al. 2008, ApJ, 684, 287

Jean-Baptiste, I., Di Matteo, P., Haywood, M., et al. 2017, A&A, 604, A106

Jofré, P., & Weiss, A. 2011, A&A, 533, A59

Koppelman, H. H., Helmi, A., & Veljanoski, J. 2018, ApJL, 860, L11

Lallement, R., Capitanio, L., Ruiz-Dern, L., et al. 2018, arXiv:1804.06060

Majewski, S. R., Nidever, D. L., Smith, V. V., et al. 2012, ApJL, 747, L37

Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2017, AJ, 154, 94

Malhan, K., Ibata, R. A., & Martin, N. F. 2018, arXiv:1804.11339

Marigo, P., Girardi, L., Bressan, A., et al. 2017, ApJ, 835, 77

Nissen, P. E., & Schuster, W. J. 2010, A&A, 511, L10

Nissen, P. E., & Schuster, W. J. 2011, A&A, 530, A15

Posti, L., Helmi, A., Veljanoski, J., & Breddels, M. 2018, A&A, 615, A70

Pouliasis, E., Di Matteo, P., & Haywood, M. 2017, A&A, 598, 66

Qu, Y., Di Matteo, P., Lehnert, M. D., & van Driel, W. 2011, A&A, 530, A10

Reid, M. J., Menten, K. M., Brunthaler, A., et al. 2014, ApJ, 783, 130

Schönrich, R., Binney, J., & Dehnen, W. 2010, MNRAS, 403, 1829

Schuster, W. J., Moreno, E., Nissen, P. E., & Pichardo, B. 2012, A&A, 538, A21

Searle, L., & Zinn, R. 1978, ApJ, 225, 357

Sommer-Larsen, J., & Zhen, C. 1990, MNRAS, 242, 10