Evidence for Two Early Accretion Events That Built the Milky Way Stellar Halo

TL;DR

This paper presents evidence for a second major accretion event, called the Sequoia, which contributed significantly to the Milky Way's stellar halo, distinct from the previously known Gaia Sausage event.

Contribution

It provides dynamical and chemical evidence for the Sequoia event as a separate accretion episode, characterizing its properties and its role in galaxy formation.

Findings

Sequoia contributed high-energy retrograde stars and globular clusters.

Sequoia's stellar mass is about 5×10^7 solar masses, with a total mass around 10^10 solar masses.

Sequoia stars have lower metallicity and distinct chemo-dynamical signatures compared to the Gaia Sausage.

Abstract

The Gaia Sausage is the major accretion event that built the stellar halo of the Milky Way galaxy. Here, we provide dynamical and chemical evidence for a second substantial accretion episode, distinct from the Gaia Sausage. The Sequoia Event provided the bulk of the high energy retrograde stars in the stellar halo, as well as the recently discovered globular cluster FSR 1758. There are up to 6 further globular clusters, including $\omega$~Centauri, as well as many of the retrograde substructures in Myeong et al. (2018), associated with the progenitor dwarf galaxy, named the Sequoia. The stellar mass in the Sequoia galaxy is $\sim 5 \times 10^{7} M_\odot$, whilst the total mass is $\sim 10^{10} M_\odot$, as judged from abundance matching or from the total sum of the globular cluster mass. Although clearly less massive than the Sausage, the Sequoia has a distinct chemo-dynamical signature. The strongly retrograde Sequoia stars have a typical eccentricity of $\sim0.6$, whereas the Sausage stars have no clear net rotation and move on predominantly radial orbits. On average, the Sequoia stars have lower metallicity by $\sim 0.3$ dex and higher abundance ratios as compared to the Sausage. We conjecture that the Sausage and the Sequoia galaxies may have been associated and accreted at a comparable epoch.

Tables

Table 1. Table 1: The kinematic, action, and orbital properties of the probable and possible GCs from Monte Carlo sampling.

Name	$(J_R,J_{\varphi },J_z)$	E	pericentre	apocentre	ecc.	incl.
	(kpc km s-1)	(km2s-2)	(kpc)	(kpc)		(deg)
FSR 1758	(${620}_{-130}^{+160},-{1250}_{-150}^{+150},{230}_{-40}^{+50}$)	$-{142500}_{-8400}^{+8100}$	${3.8}_{-0.5}^{+0.5}$	${16.2}_{-2.6}^{+3.0}$	${0.62}_{-0.01}^{+0.02}$	${146.5}_{-1.4}^{+1.4}$
NGC 3201	(${900}_{-120}^{+140},-{2800}_{-90}^{+80},{310}_{-20}^{+20}$)	$-{112300}_{-3000}^{+3200}$	${8.4}_{-0.1}^{+0.1}$	${29.3}_{-2.1}^{+2.6}$	${0.55}_{-0.02}^{+0.03}$	${152.6}_{-0.3}^{+0.3}$
$\omega$ Centauri	(${270}_{-30}^{+30},-{520}_{-30}^{+30},{100}_{-20}^{+20}$)	$-{185000}_{-800}^{+900}$	${1.5}_{-0.1}^{+0.1}$	${7.2}_{-0.2}^{+0.1}$	${0.65}_{-0.03}^{+0.02}$	${139.7}_{-2.1}^{+2.2}$
NGC 6101	(${1400}_{-210}^{+240},-{3210}_{-210}^{+200},{800}_{-50}^{+60}$)	$-{97200}_{-4400}^{+4400}$	${10.7}_{-0.5}^{+0.5}$	${41.4}_{-4.2}^{+4.7}$	${0.59}_{-0.02}^{+0.02}$	${143.1}_{-0.5}^{+0.5}$
NGC 6535	(${140}_{-30}^{+30},-{350}_{-30}^{+30},{66}_{-7}^{+8}$)	$-{207600}_{-700}^{+800}$	${1.3}_{-0.1}^{+0.1}$	${4.5}_{-0.1}^{+0.2}$	${0.56}_{-0.04}^{+0.05}$	${159.0}_{-2.4}^{+2.2}$
NGC 6388	(${80}_{-5}^{+10},-{250}_{-30}^{+40},{71}_{-1.3}^{+1.4}$)	$-{222300}_{-3300}^{+3800}$	${0.9}_{-0.1}^{+0.1}$	${3.5}_{-0.5}^{+0.2}$	${0.59}_{-0.02}^{+0.03}$	${149.2}_{-4.3}^{+3.2}$
NGC 6401	(${54}_{-7}^{+10},-{590}_{-110}^{+110},{150}_{-30}^{+30}$)	$-{194000}_{-8900}^{+8200}$	${2.5}_{-0.5}^{+0.5}$	${4.9}_{-0.8}^{+0.9}$	${0.34}_{-0.02}^{+0.02}$	${142.5}_{-1.1}^{+0.9}$

Equations

$$G$$

$$G_{\mathrm{BP}}-G_{\mathrm{RP}}$$

$${\mathcal{L}}_i={\mathcal{L}}_i^{\mathrm{clust}}+{\mathcal{L}}_i^{\mathrm{field}}.$$

$${\mathcal{L}}_i^{\mathrm{clust}}\equiv {\Sigma }^{\mathrm{clust}}(R_i)A_i^{\mathrm{clust}},{\mathcal{L}}_i^{\mathrm{field}}\equiv {\Sigma }^{\mathrm{field}}{A}_i^{\mathrm{field}},$$

$$A_i^{\mathrm{clust}}$$

$${\mathsf{S}}_i$$

$$\ln \mathcal{L}\equiv \sum _{i=1}^{N^{\mathrm{total}}}\ln {\mathcal{L}}_i,$$

$$P_i^{\mathrm{clust}}=\frac{{\mathcal{L}}_i^{\mathrm{clust}}}{{\mathcal{L}}_i^{\mathrm{clust}}+{\mathcal{L}}_i^{\mathrm{field}}}.$$

Full text

††pubyear: 2019††pagerange: Evidence for Two Early Accretion Events That Built the Milky Way Stellar Halo–Evidence for Two Early Accretion Events That Built the Milky Way Stellar Halo

1 Introduction

The Gaia Data Release 2 (DR2, Gaia Collaboration et al., 2018a) is proving transformational in the identification of substructure in the Milky Way galaxy. This is because substructure retains coherence in phase space over very long timescales (Johnston et al., 1996; Tremaine, 1999). The acquisition of kinematic data, particularly accurate stellar proper motions courtesy of the Gaia satellite, is therefore the key to unlocking the accretion history of the stellar halo. The long-term goal of understanding the building blocks of at least the stellar halo, and perhaps even the entire Galaxy, seems to be within our grasp.

Already, Gaia has provided compelling evidence for the nearly head-on collision of a Magellanic-sized dwarf galaxy with the nascent Milky Way some 8 to 10 billion years ago, the so-called ‘Gaia Sausage’ (see e.g. Belokurov et al., 2018; Myeong et al., 2018c, d; Haywood et al., 2018; Fattahi et al., 2019). This name describes the elongated shape of the structure in velocity space. The radial velocity dispersion of the Sausage stars is $\approx 180$ kms*-1*, while the azimuthal and longitudinal dispersions are only $\approx 60$ kms*-1* (see e.g., Figure 4 of Belokurov et al., 2018). The name therefore follows the long-standing scientific practice of being descriptive and informative111EV is unhappy with the name and is looking for a better one..

The aftermath of this accretion event is detectable in the inner stellar halo of the Galaxy as a giant cloud of relatively metal-rich ([Fe/H]$\gtrsim-1.5$) stars on highly radial orbits. Originally traced with nearby Main Sequence stars (Belokurov et al., 2018; Myeong et al., 2018c), the Sausage debris has now been found over a large distance range with a number of distinct tracers, including Blue Horizontal Branch stars (Deason et al., 2018; Lancaster et al., 2019) and RR Lyrae (Simion et al., 2019; Iorio & Belokurov, 2019). The characteristic property of the residue of this collision is that the orbits are eccentric with little or no net angular momentum. The debris of this event does not provide any strongly prograde or retrograde material, as befits an almost head-on collision.

Alternatively, it was proposed that an ancient major merger – dubbed ‘Gaia-Enceladus’ – could have given rise to the bulk of the retrograde stars in the halo, as well as some of the low-angular momentum debris (Helmi et al., 2018). The fundamental difference between the two hypotheses is that the ‘Gaia-Enceladus’ encompasses not just the highly eccentric component of the halo, but also the strongly retrogarde component. For example, in Helmi et al. (2018), the ‘Gaia-Enceladus’ stars have angular momentum component satisfying $-15000<J_{\phi}<150$ kms*-1* kpc*-1*, independent to the total energy, and so span a range from mild prograde through highly eccentric to strongly retrograde. The question of whether a single collision could produce such a spray of debris with different kinematical properties remains open.

In fact, the suggestion that the retrograde component of the halo may have been accreted already predates the arrival of the Gaia data by many years (see e.g., Norris & Ryan, 1989; Carollo et al., 2007; Beers et al., 2012; Majewski et al., 2012). The retrograde and peculiar globular cluster $\omega$ Centauri has also long been suspected of playing a role in the supply of retrograde stars, as it may be the stripped nucleus of a dwarf galaxy (Bekki & Freeman, 2003). The Gaia data releases have provided new samples of the retrograde halo component, which have been scoured for evidence of multiple minor mergers and accretion events (e.g., Helmi et al., 2017; Myeong et al., 2018c, b). The question therefore at issue is: did one merger event provide both the eccentric and retrograde components of the stellar halo (as in the ‘Gaia-Enceladus’ theory) or does the retrograde component have a different origin from the eccentric component (the ‘Gaia-Sausage’ theory)?

We provide a possible answer to this question in our paper, but our line of reasoning begins in a roundabout way with another unusual retrograde object, FSR 1758. This was originally discovered by Froebrich et al. (2007) as a claimed open cluster and later identified as a globular cluster (Cantat-Gaudin et al., 2018). Barbá et al. (2019) recently reported the first estimate of its distance and noticed its unusual size, using a combination of data from the DECam Plane Survey (DECaPS, Schlafly et al. 2018) and the VISTA Variables in the Via Lactea (VVV) Extended Survey, complemented with Gaia DR2. It is an extended agglomeration of stars, located at ($\ell=349^{\circ},b=3^{\circ}$) and with a heliocentric distance of $10-12$ kpc. They determined the core radius of FSR 1758 to be $\approx 10$ pc, and estimated the tidal radius to be $\approx 150$ pc. Considering its unusual size, Barbá et al. (2019) questioned whether FSR 1758 is the remnant of a dwarf galaxy or an unusually large globular cluster. Subsequently, Simpson (2019) found 3 stars in the centre of FSR 1758 with line-of-sight velocities from the Gaia Radial Velocity Spectrograph (RVS) and argued that the object has a line-of-sight velocity of $227\pm 1$ kms*-1*. Although the 3 stars are insufficient to come to a definite conclusion as regards FSR 1758’s internal velocity dispersion, nonetheless Simpson (2019) argued on the basis of its highly retrograde orbit that it is an accreted halo globular cluster.

In fact, the globular cluster (GC) datasets have been scrutinised carefully for evidence of accretion events in recent years. They have proved surprisingly powerful tracers of the merger events that build the stellar halo of the Galaxy. This was made clear in Forbes & Bridges (2010), who showed that the globular clusters that formerly belonged to the Sagittarius galaxy follow a different track in the plane of age-metallicity as compared to the bulk of the primordial or in situ clusters. Subsequently, Myeong et al. (2018d) identified a sample of 10 high eccentricity, high-energy, old halo GCs strongly clumped in action space that belonged to the ‘Gaia Sausage’ event. More speculatively, Kruijssen et al. (2019) used the age–metallicity distribution of Galactic globular clusters to reconstruct the entire merger history, claiming three substantial events. If FSR 1758 is indeed an accreted GC, then this suggests that a systematic search for companion GCs accreted in the same event, as well as other stellar debris such as substructures and tidal tails, may provide a picture of the progenitor.

In Section 2, we use Gaia’s kinematic data in combination with photometry from DECaPS to determine structural parameters and the proper motion dispersion profile for FSR 1758, and show that its declining fall-off is characteristic of a GC. The nature of FSR 1758 having been established, we search for companion GCs and stellar substructures moving on orbits of similar eccentricity and inclination that may have joined the Milky Way in the same accretion event in Section 3. Barbá et al. (2019) introduced the picturesque term Sequoia to describe the size of FSR 1758. We retain the term and slightly adapt it for our own use. In our picture, FSR 1758 is one of about five or more GCs that populated the Sequoia dwarf galaxy, whose existence was already conjectured from our stellar substructure searches (Myeong et al., 2018c, b). Its disruption brought these GCs into the Milky Way on similar orbits, as well as abundant retrograde high energy stellar substructure. We argue that the remnants of the Sequoia galaxy are dynamically distinct from the Gaia Sausage as they are retrograde, whereas the Sausage was an almost head-on collision. The dual pattern of these accretion events is evident in energy and actions, and is also shown clearly when the chemical evidence is analysed. The stars and substructures associated with the Sequoia have different mean metallicities and different abundance ratios. In Section 4, we provide estimates of the age and mass of the Sequoia galaxy and compare with the Gaia Sausage. Finally, we summarise our results in our concluding Section 5.

2 The Nature of FSR 1758

2.1 Data

First, we cross-match the positions of stars between Gaia and DECaPS, using a search radius of 0.5”. DECaPS provides roughly two magnitudes deeper photometry than Gaia, but it is saturated for bright stars and has patchy spatial coverage, so we study the union of the two datasets. We use the following combinations of $r$ and $i$ DECaPS photometric bands as a proxy for Gaia $G$ and $G_{\mathrm{BP}}-G_{\mathrm{RP}}$ (derived by comparing the magnitudes of cross-matched stars):

[TABLE]

Fig. 1 shows the surface density profiles of stars in different ranges of magnitudes, and their distributions by magnitudes at different spatial locations. By comparing the density of stars in Gaia and DECaPS datasets, we conclude that the former is reasonably complete up to $G\lesssim 20$. The subset of stars with astrometric measurements also extends roughly to $G=20$, but is less complete in the central area. We wish to include the fainter stars without astrometry in order to mitigate the bias in representation of the spatial density profile of the cluster. We chose to use stars up to $G=21$, of which roughly 70% are present in the Gaia dataset, and only 40% have astrometric measurements. During fitting, we also infer the total mass from the intrinsic (error-deconvolved) proper motion (PM) dispersion. We only use a high-quality subset of stars (marked as “good astrometry” in the figure) for this inference, ignoring all sources with astrometric_excess_noise${}>1$ or phot_bp_rp_excess_factor${}>1.3+0.06(G_{\mathrm{BP}}-G_{\mathrm{RP}})^{2}$, as suggested by Gaia Collaboration et al. (2018a).

2.2 Dynamical Modelling

We use a probabilistic model, in which the stars are drawn from a mixture of two populations: the cluster and the field. The distribution of field stars is assumed to be spatially uniform, and described by a sum of two bivariate Gaussians in the PM plane. We assume that the density of cluster stars follows a generalized King profile, also known as the limepy family of models (Gieles & Zocchi, 2015), which has been shown to adequately describe realistic globular clusters (Hénault-Brunet et al., 2019). It has the following free parameters: mass $M$, scale (core) radius $R_{c}$, dimensionless potential depth at the centre (King parameter) $W_{0}$, and truncation parameter $g$ controlling the density profile in the outer parts. For the models in this section, we do not assume any particular relation between the total cluster mass and the number of observed cluster members $N^{\mathrm{clust}}$ (this relation is examined in the next section). Rather, the total mass of the cluster manifests itself only kinematically, through the overall amplitude of velocity dispersion. We assume a Gaussian distribution for the PM of cluster stars, centred around its mean PM, and with a spatially-variable width.

We measure the parallax distribution of field stars directly from the data using stars outside the central 5 arcmin and represent it by a mixture of three Gaussian components, with parameters fixed throughout the rest of the modelling. The intrinsic (error-free) parallaxes of cluster stars are assumed to be equal to the inverse distance to the cluster (fixed to $D=10$ kpc), plus the constant zero-point parallax offset $-0.03$ mas (Gaia Collaboration et al., 2018a). We do not use any information about colours and magnitudes for membership determination, because they are severely affected by spatially-variable reddening.

The likelihood of observing a star $i$ at a given distance from the cluster centre, with or without further astrometric information, is given by a sum of likelihoods of the two alternative hypotheses:

[TABLE]

Let $N^{\mathrm{clust}}$ be the (unknown) total number of cluster members in our sample of stars within the radius $R_{\mathrm{max}}$. We choose $R_{\mathrm{max}}\gg R_{c}$, so that all possible cluster members are included, and normalize the cluster surface density $\Sigma^{\mathrm{clust}}(R)$ so that $\int_{0}^{R_{\mathrm{max}}}\Sigma^{\mathrm{clust}}(R)\;2\pi\,R\,\mathrm{d}R=N^{\mathrm{clust}}$. The remaining $N^{\mathrm{field}}\equiv N^{\mathrm{total}}-N^{\mathrm{clust}}$ observed stars are then attributed to the field population with a spatially-uniform density $\Sigma^{\mathrm{field}}\equiv N^{\mathrm{field}}/(\pi\,R_{\mathrm{max}}^{2})$. Then the likelihoods of a given star to belong to either population are

[TABLE]

where the factors $A^{\mathrm{clust}}_{i},\>A^{\mathrm{field}}_{i}$ are unity for stars without astrometric measurements, or describe the likelihood of measuring the observed parallax and PM, given the intrinsic distribution functions of either population, convolved with measurement uncertainties. For the cluster population, this factor is a product of the parallax likelihood and the PM likelihood (ignoring correlations between them):

[TABLE]

where $\mathcal{N}$ is the uni- or bivariate normal distribution, $\varpi$ is the parallax, $\boldsymbol{\mu}\equiv\{\mu_{\alpha},\mu_{\delta}\}$ is the PM with associated measurement uncertainties $\epsilon$, $\rho_{i}$ is the correlation coefficient between the two components of PM uncertainty matrix, $R_{i}$ is the distance of the star from the cluster centre, $\sigma_{\mu}(R)$ is the spatially-dependent intrinsic PM dispersion of cluster stars, whose amplitude is proportional to the square root of the cluster mass. We use only a subset of stars with reliable PM measurements outside the central 2 arcmin for inferring the intrinsic dispersion, since the stars in the centre may be affected by crowding. The PM of remaining (mostly faint) stars are still used to determine the membership, but in doing so, we use a conservative value $\sigma_{\mu}=0$. For the field population, the expressions are similar, but involve several Gaussian components (two for PM and three for parallax), with the intrinsic PM dispersion being a spatially-constant symmetric $2\times 2$ matrix rather than a single spatially-varying quantity.

The fitting procedure optimizes the model parameters ($N^{\mathrm{clust}},\boldsymbol{\mu}^{\mathrm{clust}},\boldsymbol{\mu}^{\mathrm{field}}$, covariance matrices of the field population, parameters of the cluster density profile) to maximize the total log-likelihood

[TABLE]

via Markov Chain Monte Carlo as implemented in the emcee package (Foreman-Mackey et al., 2013). We then evaluate the posterior probability of membership for each star as

[TABLE]

The total number of cluster members is $N^{\mathrm{clust}}=\sum_{i}P^{\mathrm{clust}}_{i}$. We stress that we do not make hard cuts in any of the observed quantities (parallax, PM, radius) to separate the cluster and the field populations. For stars with astrometric measurements, the membership probability distribution is strongly bimodal (the two populations are well separated), while for stars without astrometry, the membership probability smoothly drops with radius from $\sim 0.5$ in the centre down to zero at large radii. In total, we have $N^{\mathrm{clust}}\simeq 7500$ member stars, of which only $\sim 1600$ are astrometrically selected, and only $\sim 350$ are used in the dynamical mass determination through the intrinsic PM dispersion.

Fig. 2 shows the inferred cluster density and the PM dispersion profiles. The best-fit parameters of the generalized King model are: King radius of 3 arcmin ($\sim 8.5$ pc for the assumed distance $D=10$ kpc), King parameter $W_{0}\simeq 3$, truncation radius $\sim 15$ arcmin, and truncation parameter $g\simeq 1$. The core radius (defined as the projected distance from the cluster centre where the surface density drops to 1/2 its central value) is $\sim 6.5$ pc, and the half-light radius (the projected distance enclosing half of the cluster stars) is $\sim 9$ pc. The core radius is somewhat smaller than determined by Barbá et al. (2019) from their fit to the DECaPS photometric sample (without considering astrometry). On the other hand, the truncation radius ($\sim 50$ pc) is three times smaller than found in that paper. We stress that the density profile in the outer parts, and in particular the truncation radius, is determined mainly by Gaia astrometry, so is more reliably constrained than just using the photometry alone.

Overall, the modelling procedure makes good use of both the deeper DECaPS photometry in the cluster centre and the Gaia astrometry in the outer parts. However, the inferred cluster mass (equivalently, the intrinsic PM dispersion) appears to be rather high, $M\simeq(7\pm 1)\times 10^{5}\,M_{\odot}$, compared to the photometric model of the next section. We stress that the width of the intrinsic PM distribution is inferred by convolving it with the measurement uncertainties and comparing the error-broadened distribution with the actually observed one. Hence, it strongly depends on the reliability of uncertainty estimates $\epsilon_{\mu}$ of stars in the Gaia dataset. Even for the high-quality subsample, these errors are in the range $0.1-0.3$ mas yr*-1*, comparable to or exceeding the inferred value of intrinsic PM dispersion. It is known that the formal uncertainties in Gaia PMs are underestimated (e.g., Gaia Collaboration et al., 2018a). We multiplied the uncertainties quoted in the catalogue by a correction factor 1.1, as suggested in that paper, before running the fit. If we instead increase the uncertainties by a factor 1.3, this reduces the PM dispersion by a third, bringing it into agreement with the total cluster mass estimated from photometry.

Finally, to check if the PM dispersion profile could be biased by the assumed parametric form of the generalized King model, we also determined it non-parametrically from the subset of high-quality stars classified as cluster members with $\geq 80\%$ probability. We used the method of Vasiliev (2018), representing $\sigma_{\mu}(R)$ as a cubic spline with the values at four control points adjustable during the fit, while taking into account spatially correlated measurement errors. The result, shown in the right panel of Fig. 2, agrees reasonably well with the parametric profile, but is higher in the very centre. This is likely caused by crowding issues, and for this reason we have excluded the stars in the central two arcmin from the high-quality sample used to determine the PM dispersion in the parametric fit. In any case, the PM dispersion profile appears to be declining with radius, which is natural to expect for a globular cluster, but not for a dwarf galaxy.

2.3 Photometric modelling

Since almost all stars with magnitudes $G\lesssim 19-20$ have astrometric measurements and are well separated into the cluster and the field populations, we may use the sample of astrometrically confirmed members to determine the total cluster mass. The idea is to compare the distribution of stars by magnitudes with those of several other clusters with similar colour-magnitude diagrams (CMDs), for which the total mass and the distance are known. By normalizing the number of stars brighter than the given absolute magnitude by the cluster mass, we constructed the cumulative magnitude distribution profiles for a dozen globular clusters in the range of metallicities $-2\lesssim\mathrm{[Fe/H]}\lesssim-1$. The masses and distances of these clusters are taken from Baumgardt et al. (2019), who used a large library of $N$-body simulations and a variety of observational constraints to measure the masses. Three clusters have very similar CMDs to FSR 1758, particularly with regards to the location of blue horizontal branch (BHB): NGC 6205 (M 13), NGC 6254 (M 10) and NGC 6656 (M 22). They all have metallicities $\simeq-1.5$ and masses $(2-5)\times 10^{5}\,M_{\odot}$.

Fig. 3, left panel, shows the composite CMDs of these three clusters, with the astrometrically detected members of FSR 1758 overplotted by black crosses. The measured colour $G_{\mathrm{BP}}-G_{\mathrm{RP}}$ and the $G$-band magnitude are de-reddened using the coefficients given in Table 2 of Gaia Collaboration et al. (2018b), and the observed magnitude is converted to the absolute magnitude. By matching the stars of FSR 1758 to the empirically determined isochrones, we infer the reddening coefficient $E(B-V)=0.73$ and the distance $D=10$ kpc, in reasonable agreement with Barbá et al. (2019). Stars in the lower part of the red giant branch (RGB) are systematically offset to the left from the isochrone, but this is expected for such a dense and highly reddened region. If we use DECaPS $r$ and $i$ bands instead of Gaia $G_{\mathrm{BP}}-G_{\mathrm{RP}}$, as per eqn 1, the scatter and offset of member stars from the isochrone curve at faint magnitudes are substantially reduced. The upper part of the RGB and the BHB of FSR 1758 match well the location of these features in the other three clusters.

The right panel of Fig. 3 shows the cumulative number of stars as a function of magnitude, normalized by the mass of each cluster. The three clusters listed above have similar profiles, and in fact almost all other clusters also follow the same trend (we show additionally the cluster NGC 6752, which has a somewhat higher metallicity). By matching, we infer the mass of FSR 1758 to be $\sim 2.5\times 10^{5}\,M_{\odot}$, with an uncertainty $\lesssim 20\%$.

As noticed by Simpson (2019), three bright stars in the centre of FSR 1758 are actually present in the Gaia RVS sample, having values of line-of-sight velocity around 227 km s*-1*. Simpson (2019) additionally reported a fourth star with a very similar line-of-sight velocity at a projected distance $\sim 0.6^{\circ}$ from the cluster. Given that its PM, $G$ magnitude and colour are all very close to those of the three other stars, it is unlikely to be a field star, but it also cannot be a current cluster member, being more than twice as far as the inferred cutoff radius. This star may have been tidally stripped from the cluster.

With only three stars, it is not possible to put strong constraints on the internal velocity dispersion $\sigma$. The measured values are very close to each other and consistent with zero intrinsic scatter, but values of $\sigma$ up to 5–10 km s*-1* are not strongly excluded (Fig 4). These values are also consistent with the photometrically estimated mass, which corresponds to the central velocity dispersion $\sim 4$ km s*-1*, and even with the higher values inferred from PM, although these are less reliable.

3 Tracers of the Sequoia Galaxy

FSR 1758 has a very retrograde, eccentric orbit (e.g., Simpson, 2019). An accretion origin for the strongly retrograde components of the stellar halo has long been suspected (e.g., Quinn & Goodman, 1986; Norris & Ryan, 1989; Carollo et al., 2007; Beers et al., 2012). Analyses of the Gaia data (Helmi et al., 2017; Myeong et al., 2018c, b) have convincingly shown that the highest energy stars in the halo are typically retrograde, reinvigorating these earlier suspicions. This rotational asymmetry could be the consequence of dynamical interaction between accreted satellites and the Milky Way. For example, Quinn & Goodman (1986) and Norris & Ryan (1989) show that retrograde accretion events – especially with some inclination – experience much less drag than prograde ones. So, the effect of many random infalls can produce an the overall rotational asymmetry at high energy. A sharp feature traced by the retrograde stellar substructure has already been seen in energy-action space with a specific metallicity range. This is especially obvious in the subpanels of metallicity $-1.9<\mathrm{[Fe/H]}<-1.5$ in Figure. 2 of Myeong et al. (2018c). This is a clear signature of an individual event of retrograde accretion in the past, distinct from other, numerous, past accretions.

Fig. 5 shows the Milky Way GCs from Vasiliev (2019), together with the substructures found by Myeong et al. (2018b), in the plot of scaled action. GCs on retrograde orbits, including FSR 1758, lie on the far left-hand side of the plot. Specifically, the horizontal axis is the (normalized) azimuthal action $J_{\phi}/J_{\mathrm{tot}}$), while the vertical axis is the (normalized) difference between the vertical and radial actions ($(J_{z}-J_{R})/J_{\mathrm{tot}}$). Colour represents the radius of the circular orbit with the same energy, and so gives an idea of the typical distances probed by an object. The observational uncertainties such as the distance, line-of-sight velocities and proper motions are Monte Carlo sampled, with error ellipse transforming in action space to distended shapes. The gravitational potential used to represent the Milky Way is the one recommended as the best amongst the suite studied by McMillan (2017). It is an axisymmetric model with bulge, thin, thick and gaseous disks and an NFW halo.

The portion of the plot occupied by FSR 1758 (coloured green) overlaps with a number of GCs, in particular $\omega$ Centauri. Using a ‘Friends-of-Friends’ clustering algorithm in this scaled action space, we identify 6 GCs that form an agglomeration. They are FSR 1758, NGC 3201, $\omega$ Centauri (NGC 5139), NGC 6101, NGC 5635, and NGC 6388. All 6 are listed in Table 1, which gives their actions, energies and orbital characteristics. Also shown on Fig. 5 are all the retrograde stellar substructures identified by Myeong et al. (2018b). These are depicted as irregularly shaped polygons that include all the stars believed to be members. Apart from Rg5 and Rg7, it is striking that all the retrograde substructures overlap with our group of GCs associated with FSR 1758. Fig. 6 shows the same objects, but now plotted in the plane of eccentricity and inclination. We see that the GCs listed in Table 1 are restricted to a narrow range of inclinations ($140^{\circ}-160^{\circ}$) and eccentricities ($e\sim 0.6$). We remark that additionally NGC 6401 may be associated with the group, at least as judged by inclination. However, its eccentricity is somewhat less than the other members. As a possible member, the orbital properties of NGC 6401 are also listed in the lower part of Table 1. It is noteworthy that the inclination range of our group appears to be distinct from the orbital plane of the Magellanic system (e.g., D’Onghia & Lake, 2009; Nichols et al., 2011) or the plane of Milky Way satellites (Kroupa et al., 2005) known to be near perpendicular to the Galactic plane.

In addition to the GCs (Vasiliev, 2019) and retrograde stellar substructures (Myeong et al., 2018b), a retrograde stellar stream, GD-1 (Grillmair & Dionatos, 2006) is also marked on Fig. 5, 6 and 7 based on a representative six-dimensional phase space information from Webb & Bovy (2019). Interestingly, the orbital inclination, (normalized) azimuthal action and energy of GD-1 appear to be comparable to our group of GCs and retrograde substructures. But, its other action components and orbital eccentricity noticeably differ from our group. According to the complex morphology of GD-1 (see e.g., Price-Whelan & Bonaca, 2018; Malhan et al., 2019), there is a possibility that the current orbital characteristics of GD-1 may not be a good reflection of its past state or its progenitor or its parent dwarf galaxy. Nonetheless, it will require more detailed investigation to find out any potential connection between our group and GD-1.

There have been long-held suspicions over the intriguing and anomalous globular cluster, $\omega$ Centauri. Normal globular clusters often have multiple populations, but show homogeneous abundances in heavy elements (such as calcium and iron) together with variations in light elements (such as the oxygen-sodium anticorrelation). However, $\omega$ Centauri hosts multiple stellar populations with different heavy element abundances enriched by supernovae, as well as spreads in the light elements (see e.g., Lee et al., 1999; Bedin et al., 2004; D’Antona et al., 2011; D’Orazi et al., 2011; Joo & Lee, 2013). This necessitates the existence of different channels for enrichment to account for the chemical peculiarities of $\omega$ Centauri (e.g., Bekki & Norris, 2006; Romano et al., 2007; Marcolini et al., 2007). Already Freeman (1993) and Bekki & Freeman (2003) suggested it is the nucleus of a stripped dwarf galaxy, inspired by its very bound retrograde orbit. A number of authors have argued that the retrograde components in the stellar halo may be related to the disruption of $\omega$ Centauri (e.g., Dinescu, 2002; Brook et al., 2003; Majewski et al., 2012; Helmi et al., 2017; Myeong et al., 2018c, b). We amplify this hypothesis here by associating it with the Sequoia Event. Either $\omega$ Centauri is the remnant core of the Sequoia galaxy, or it was the largest GC member. For dwarf galaxies in cored haloes, the nucleus may be completely dissolved by the merging process, leaving only the GCs and stellar debris. For nucleated dwarfs or dwarfs in cusped haloes, the nucleus can survive intact, even if the outer parts are stripped. Whichever picture is correct, it is reasonable to conclude that the Sequoia dwarf once possessed an entourage of up to 7 GCs (5 probable and 2 possible).

Fig. 7 shows the distribution in energy and normalised azimuthal action. We show the azimuthal action of $\omega$ Centauri with its $3\sigma$ uncertainty as dashed vertical lines. If inclination is roughly preserved under the action of dynamical friction for these strongly retrograde mergers, then the dashed lines will delineate GCs and substructure associated with the $\omega$ Centauri, and hence the Sequoia. We see that FSR 1758, as well as the objects in Table 1, all lie within this band. GCs with higher energy were stripped earlier and/or composed the trailing tail. In this picture, FSR 1758 was one of the most massive GCs in the precursor dwarf galaxy, residing near the centre of the progenitor. So, it has ended up with comparable, if higher, energy. Objects with lower energy may have come from the disruption of the leading tail. This includes NGC 6535 (and possibly NGC 6388), which have been deposited closer to the centre of the Milky Way on tighter orbits.

From energy arguments, the retrograde substructures (Myeong et al., 2018b) are more likely to be the remnant debris of the Sequoia than direct tidal debris from $\omega$ Centauri. The retrograde stellar substructures are stripped well before the progenitor became severely destroyed. In this case, if its core is indeed $\omega$ Centauri, we do not expect the retrograde substructures to show the unique chemical abundance patterns observed from $\omega$ Centauri (e.g., Na–O and Mg–Al patterns, Ba overabundance; cf. Navarrete et al., 2015). These may have been imprinted on $\omega$ Centauri by subsequent accretion of gas onto the nucleus, after the stripping process removed the retrograde substructures. It would be interesting to revisit Kapteyn’s Moving Group (Wylie-de Boer et al., 2010) or the $\omega$ Centauri Moving group (Meza et al., 2005) – which were disproved to be $\omega$ Centauri’s tidal debris based on the chemical abundances (Navarrete et al., 2015) – for the possibility that they are remnant debris of the Sequoia. By contrast, the freshly discovered Fimbulthul stream (Ibata et al., 2019) is probably stripped from $\omega$ Centauri itself rather than Sequoia and in this case should show the unique abundance patterns.

The age–metallicity relation of the GCs is shown in Fig. 8 based on data complied by Kruijssen et al. (2019). Among five member GCs estimates four GCs (except NGC 6388) form a distinct track in the age–metallicity relation that is different from the Milky Way’s in-situ GCs. This is most visible in the upper panel of Fig. 8. This branching is similar to what has already been seen for GCs associated with major accretion events – in particular, the Sagittarius (Sgr) GCs (Terzan 7, Terzan 8, Arp 2, Pal 12, NGC 4147222As will be discussed later, this cluster is unlikely to belong to the Sgr group, based on its kinematics., NGC 6715, and Whiting 1, Forbes & Bridges, 2010, and marked as diamonds) or the Sausage GCs (see e.g., Myeong et al., 2018d, and marked as downward-pointing triangles). The track of the Sausage GCs stands out clearly especially at the lower panel of Fig. 8. This track in the age–metallicity plot of Fig. 8 is remarkable especially since the original membership of the GCs is established purely based on dynamical information. This is independent evidence of the extragalactic origin of the member GCs from an individual merger event.

Halo stars can be identified effectively based on chemical abundances such as [Al/Fe] and [Mg/Fe] (Hawkins et al., 2015). With APOGEE DR14 (Abolfathi et al., 2018), Mackereth et al. (2019) showed that halo stars with high eccentricity orbits tend to have lower [Mg/Fe] on average compared to the rest of the halo stars. In Fig. 9, we show more specifically that the highly radial Gaia Sausage remnant stars ($e\sim 0.9$ with $\lvert J_{\phi}/J_{\mathrm{tot}}\rvert<0.07$ and $(J_{z}-J_{R})/J_{\mathrm{tot}}<-0.3$) and the high energy retrograde substructure stars ($e\sim 0.6$ with $J_{\phi}/J_{\mathrm{tot}}<-0.5$ and $(J_{z}-J_{R})/J_{\mathrm{tot}}<0.1$) have clearly different [Fe/H] distribution and different abundance patterns. The Sausage remnants show a metallicity distribution function peak at [Fe/H]$=-1.3$, whereas the high energy retrograde stars are more metal poor, with the peak at [Fe/H]$=-1.6$ (both in good agreement with Myeong et al., 2018c; Matsuno et al., 2019). While the metallicitiy distributions of the Sausage and Sequoia stars overlap, at fixed [Fe/H], the two galaxies show distinct patters in the abundance of alpha elements. For example, at [Fe/H]$\sim-1.5$, the Sequoia debris are clearly more enhanced in Al compared to the Sausage. Such differences in the abundances provide additional evidence that the accretion event that made the high energy retrograde stars in the halo is also chemically different from the Sausage event. Interestingly, their chemical characteristics – the Sausage having higher [Fe/H] and lower abundance ratios compared to the Sequoia progenitor – are in line with the trend observed by Mackereth et al. (2019) from the EAGLE simulations (Schaye et al., 2015).

Matsuno et al. (2019) searched through the SAGA database, which contains $\sim 880$ metal-poor stars with [Fe/H] $<-0.7$. They also found that the high energy retrograde stars are clearly distinct from the stars of the Sausage which dominates the inner halo (Belokurov et al., 2018). They reported that the ‘knee’ in the abundance and metallicity plane differs by about 0.5 dex (at [Fe/H]$\sim-2$ for the Sausage and $\sim-2.5$ for the retrograde stars) which is another indication of their different origin.

Our hypothesis is distinct from the proposal of Helmi et al. (2018), who made a broad selection based on the azimuthal action $-1500$ kpc km s*-1* ${}<J_{\phi}<150$ kpc km s*-1* only, as opposed to using additional integrals of motion. Helmi et al. (2018) used this sample to suggest the ‘Gaia-Enceladus’ accretion event. In our picture, this sample contains stars belonging to both the Gaia Sausage and the Sequoia. Our hypothesis is more closely related to the work of Mackereth et al. (2019), who divided halo stars according to eccentricity and showed that the low and high eccentricity groups have different abundance ratios and probably different origin. For the low eccentricity group, Mackereth et al. (2019) suggested they are likely to be a mixture of in situ halo stars and many smaller accreted materials – which includes the Sequoia debris. Matsuno et al. (2019) finding of a different ‘knee’ is important corroboratory evidence of this as well.

4 The Sequoia and the Sausage

The progenitor of the Sausage had a total mass in stars and dark matter \gtrsim 10^{11}$$M_{\odot} (Belokurov et al., 2018). This was also derived in Fattahi et al. (2019) using cosmological hydrodynamical simulations (Grand et al., 2017), and is consistent with the estimates from Mackereth et al. (2019) and Vincenzo et al. (2019). Myeong et al. (2018d) identified 10 GCs associated with the Sausage event from one of the earliest kinematics dataset of 75 Milky Way GCs based on the Gaia DR2 (Gaia Collaboration et al., 2018c). This Sausage GC membership has been revised with the more complete catalogue of 150 Milky Way GCs (Vasiliev, 2019) and also with the age–metallicity relation (see e.g., lower panel of Fig. 8). The original identification by Myeong et al. (2018d) was also limited to GCs with noticeably high energy only. These authors used the total energy of the young halo (YH) clusters (see e.g., Zinn, 1993; Mackey & Gilmore, 2004; Mackey & van den Bergh, 2005, for more details on the GC classification) known at that time as a reference energy in order to concentrate on the old halo (OH) GCs with the most certain ex situ origin, judging from their noticeably high total energy. As we now have better knowledge on the typical chemo-dynamical characteristics of the Sausage GCs (e.g., $e\sim 0.9$, $(J_{z}-J_{R})/J_{tot}\sim-0.6$ with age–metallicity branching from Fig. 5 of Vasiliev (2019) and Fig. 8 in this work), we can revise the search in the larger catalogue of Milky Way GCs, without any artificial energy cut. We note that only the GCs with full kinematics, age and metallicity information have been examined. This gives up to 21 potential GCs showing typical chemo-dynamical characteristics of the Sausage, specifically NGC 362 333YH clusters (see e.g., Zinn, 1993; Mackey & Gilmore, 2004; Mackey & van den Bergh, 2005), 1261 ${}^{\ref{ftn:note}}$, 1851, 1904, 2298, 2808, 4147 ${}^{\ref{ftn:note}}$444Note that NGC 4147 is among the potential Sausage GCs. NGC 4147 has previously been suggested as a Sgr GC (e.g., Bellazzini et al., 2003; Forbes & Bridges, 2010), while there are studies suggesting no connection with Sgr dwarf (e.g., Law & Majewski, 2010). Here, we consider NGC 4147 to be a potential Sausage GC as the orbital characteristics of NGC 4147 are very different from the Sgr dwarf and other Sgr GCs, while they are similar to other Sausage GCs., 4833, 5286, 5694, 6544, 6584 ${}^{\ref{ftn:note}}$, 6712, 6779, 6864, 6934 ${}^{\ref{ftn:note}}$, 6981 ${}^{\ref{ftn:note}}$, 7006 ${}^{\ref{ftn:note}}$, 7089, Pal 14 ${}^{\ref{ftn:note}}$, Pal 15. Note that there are 8 YH classified GCs among the set with red horizontal branch (HB) morphology (see e.g., Zinn, 1993; Mackey & Gilmore, 2004; Mackey & van den Bergh, 2005). Since the Sausage progenitor was considerably massive, its original GCs might have HB morphology similar to the Milky Way OH clusters (see Mackey & Gilmore, 2004; Mackey & van den Bergh, 2005, for more detail), while some of the younger GCs with red HB morphology might have been formed during the wet merger of the Sausage and the Milky Way (see e.g., Renaud, 2018), or have been acquired by the Sausage from separate accretion events before it merged with the Milky Way. Interestingly, this number, 13, of Sausage GCs with OH classification is in agreement with the suggested number of expected GCs originally formed in a major merger satellite of the Milky Way (see e.g., Kruijssen et al., 2019).

The Sequoia event has at least 4 GCs, excluding $\omega$ Centauri (as is right if it is the stripped core). The total stellar mass of GCs is known to be correlated with the (halo) mass of the host, albeit with some scatter. Based on the current and the initial masses of the GCs from Baumgardt et al. (2019) and using the fractional mass $\langle\eta\rangle=M_{\mathrm{GCs}}/M_{\mathrm{halo}}\sim 4\times 10^{-5}$ derived by Hudson et al. (2014), we can estimate the progenitor mass. In the case of the Gaia Sausage with 13 OH clusters, the progenitor mass is at least 1\times 10^{11}$$M_{\odot} (from the current mass of the GCs) or up to 4\times 10^{11}$$M_{\odot} (from the initial mass of the GCs). Hudson et al. (2014) also provides a fractional value for the GC mass in terms of the stellar halo mass, which yields corresponding lower bounds on the stellar mass of the Sausage progenitor of 5\times 10^{8}$$M_{\odot} or up to 5\times 10^{9}$$M_{\odot} respectively. The abundance matching relation of Garrison-Kimmel et al. (2014) gives similar estimates as well. Since our estimates are not redshift corrected, we note that they are overestimates of the actual progenitor mass at the merger time.

In comparison, the mass of the Sequoia galaxy is at least 1\times 10^{10}$$M_{\odot} (from the current mass of the GCs) or 5\times 10^{10}$$M_{\odot} (from the initial mass555For the mass fraction $M_{\mathrm{current}}/M_{\mathrm{initial}}$ for FSR 1758, we used the mean mass fraction of the other member GCs. of the GCs). Valcarce & Catelan (2011) suggested the mass of the progenitor of $\omega$ Centauri may be as high as 10^{10}$$M_{\odot} from chemical evolution modelling of its multiple populations which is in broad agreement. The lower bound of the stellar mass from abundance matching gives 5\times 10^{6}$$M_{\odot} (current mass of GCs) or 7\times 10^{7}$$M_{\odot} (initial mass), while the relation from Hudson et al. (2014) gives higher stellar mass estimates by a factor of two. The mass-metallicity relation of Kirby et al. (2013) with a metallicity of $-1.6$ gives a broadly consistent stellar mass of 2\times 10^{7}$$M_{\odot}, but if we take account of the redshift evolution of the relation, this could be larger (see e.g., Ma et al., 2016). For example, for $z=1.3$, the relation gives 1.7\times 10^{8}$$M_{\odot} which is again comparable to other estimates.

Although less massive then the Gaia Sausage, the Sequoia was a notable accretion in the evolutionary history of the Milky Way. In terms of the mass, the Fornax dwarf spheroidal could be a rough representation of the Sequoia progenitor. The fact that the Fornax dSph hosts a comparable number of GCs is also a similarity. Among its six GCs, Fornax 4 has been considered for possible ex situ origin (see e.g., Buonanno et al., 1999; van den Bergh, 2000), leaving Fornax with five probable in situ GCs. This is comparable to our estimates of the current number of member GCs of the Sequoia galaxy.

It is noteworthy that NGC 3201 is one of the probable members of the Sequoia event. It has already been pointed out that NGC 3201 is potentially associated with the S1 stellar stream (Myeong et al., 2018a; O’Hare et al., 2018). In fact, NGC 3201 and the S1 stream have almost identical actions and energy. Myeong et al. (2018a) identified the S1 stream as a stellar remnant of an accreted dwarf, and inferred its progenitor mass to be in order of \sim 10^{10}$$M_{\odot} based on the library of accretion events using minor merger N-body simulations (Amorisco, 2017). Interestingly, in this study, we already derived the lower bound mass of the Sequoia to be in order of \sim 10^{10}$$M_{\odot} as well, based on completely different methods. Such agreement from independent approaches adds credence to the authenticity of this likely the same event discovered from two separate studies.

Now, the infall time for the S1 stream’s progenitor is $\gtrsim 9$ Gyr (Myeong et al., 2018a). If the S1 progenitor is indeed the Sequoia, this provides a consistent picture, as the youngest GC associated with the event (NGC 3201) is $\sim 11$ Gyr. In fact, the S1 stream and NGC 3201 bracket the range of possible infall times between 9 and 11 Gyr. Interestingly, Marcolini et al. (2007) suggested that $\omega$ Centauri is a remnant of a dwarf spheroidal galaxy accreted $\sim 10$ Gyr ago based on their study with hydrodynamical and chemical modelling. This is in a good agreement with our range. Also, we note that this suggested age range (9 Gyr to 11 Gyr) agrees very well with the suggested infall time of the Gaia Sausage itself. There is a possibility that the Sausage and the Sequoia galaxies were accreted at a comparable epoch. This suggests that here may have been a global association between them – perhaps the Sequoia was a satellite galaxy of the Sausage? The currently observed higher angular momentum of the Sequoia debris could reflect its binary orbital velocity at the time of accretion. More generally, any small differences in the initial condition at the early stage of the infall can easily cause the current difference between the orbital characteristics of the Gaia Sausage (highly radial with very low azimuthal action) and the Sequoia event (clearly retrograde).

5 Conclusions

The starting point of our investigation is an unusual object FSR 1758 (Froebrich et al., 2007). Its enigmatic nature was recently pointed out by  Barbá et al. (2019), who raised the question as to whether FSR 1758 is an unusually large globular cluster or a dwarf galaxy remnant. Using Gaia data, we derived its proper motion dispersion profile, which is strongly declining, and so we concur with Simpson (2019) that FSR 1758 is an accreted, retrograde globular cluster. Our modelling suggests that FSR 1758 has a half-light radius of $\sim 9$ pc and a baryonic mass of $\sim 2.5\times 10^{5}\,M_{\odot}$ with an uncertainty $\lesssim 20\%$.

It is natural to look for other retrograde globular clusters with similar actions as FSR 1758, which may have fallen in to the Milky Way at the same merger. This led us to the identification of Sequoia Event, which was already conjectured from our studies of stellar substructures (e.g., Myeong et al., 2018c, b). Other investigators before us (e.g., Quinn & Goodman, 1986; Norris & Ryan, 1989; Carollo et al., 2007; Beers et al., 2012; Majewski et al., 2012) have concurred that the highly retrograde parts of the stellar halo are most likely accreted.

The Sequoia Event is distinct from other known accretions, particularly the Gaia Sausage (Belokurov et al., 2018; Myeong et al., 2018c, d). It has been seen in three tracers. First, there are at least 6 globular clusters, packed in action space around FSR 1758. These include $\omega$ Centauri itself, the most massive of the Milky Way globular clusters. Based on its kinematics and its spread of stellar ages, metallicities and abundances, this has long been suggested as the stripped core of a dwarf galaxy remnant (e.g., Bekki & Freeman, 2003; D’Antona et al., 2011). It may therefore be the remnant of the Sequoia galaxy. However, it is also possible that $\omega$ Centauri may have been a core globular cluster of a now wholly destroyed progenitor. Whichever hypothesis is correct, it still remains the case that FSR 1758 is also one of the Sequoia’s largest globular clusters. Of the agglomeration of 6 globular clusters, 4 have existing age and metallicity estimates (see e.g., Kruijssen et al., 2019, and the references therein). They form a distinct track in the age–metallicity relation, different from the Milky Way’s in situ globular clusters. Additionally, this track shows evidence of an offset from the track of the Sausage globular clusters (Myeong et al., 2018a), supporting the identification of a separate event.

Secondly, the Sequoia Event is discernible in the retrograde stellar substructures, which are clearly visible in energy and action space (Myeong et al., 2018b). They form a separate grouping from the bulk of the Gaia Sausage, which has close to zero net angular momentum. This is clearly shown in Fig. 10, where the morphology of the contours in the high energy region shows the pattern of bimodal accretion tracks. The distribution of the azimuthal action for the stars with high energy (e.g., $E>-1.1\times 10^{5}$ km2/s2) shows the existence of this extra retrograde component, clearly separated from the Sausage at zero angular momentum. Myeong et al. (2018c, b) showed the signal of this extra component is concentrated at a specific range of metallicity ([Fe/H] $\sim-1.6$). Thus, the Sequoia Event is also distinct from the Gaia-Enceladus structure (Helmi et al., 2018), which appears to combine parts of the Gaia Sausage and the Sequoia.

Thirdly, the very metal poor stars that are retrograde also have a chemical signature in the abundance and metallicity plane that is distinct from both the Sausage and the overall halo. This was already hinted at in Mackereth et al. (2019), who used APOGEE data release 14 (DR14, Abolfathi et al., 2018) to demonstrate that the retrograde halo stars have lower [Mg/Fe] compared to the rest of the halo. The argument was further substantiated by Matsuno et al. (2019), who found evidence that the knee in the abundance and metallicity plane occurs at different spots separated by 0.5 dex for the stars in the Sausage and in the retrograde component or Sequoia. We have provided further evidence from APOGEE DR14 that the metallicity distribution and the abundance patters of Sausage and Sequoia stars are different (see Fig. 9). The peak of the metallicity distribution function of the Sausage is higher with [Fe/H]${}=-1.3$ as compared to the Sequoia at [Fe/H]${}=-1.6$. The abundance ratios of Sausage stars are lower than Sequoia stars.

These three lines of evidence all argue for two different accretion events. A number of arguments (abundance matching, mass in globular clusters today) suggest that the Sausage progenitor had a total mass of $\sim 1-5\times 10^{11}$ $M_{\odot}$, whilst the Sequoia had a mass of $\sim 1-5\times 10^{10}$ $M_{\odot}$. In terms of stellar mass, the Sausage weighs in at $\sim 5-50\times 10^{8}$ $M_{\odot}$, whilst the Sequoia is $\sim 5-70\times 10^{6}$ $M_{\odot}$. Modern day analogues would be the Large Magellanic Cloud and Fornax dSph, respectively. The infall time is somewhat comparable, and so the two progenitors may have been a binary pair or association.

Different tracers (e.g., globular clusters, retrograde stellar substructures) stripped from the progenitor at different times now occupy different portions of energy and action space. They are like stepping stones that enable us to recreate the history of the event and trace out the time evolution of the disruption of the progenitor. To carry out such a re-creation, we need to be certain which substructures can be definitely associated with the Sequoia event. Here, detailed studies of abundances of stars with high resolution spectroscopy can play a crucial role. The chemical signature of the Sequoia event is seemingly evident both in Fig. 9, as well as in Mackereth et al. (2019) and Matsuno et al. (2019). The current limitation is the size or the quality of the sample and so dedicated medium or high resolution spectroscopic follow-up study is essential.

acknowledgements

We thank the referee for valuable comments and suggestions, which helped to improve the quality of the publication. GCM thanks the Boustany Foundation, Cambridge Commonwealth, European & International Trust and Isaac Newton Studentship for their support of his work. EV acknowledges support from the Science and Technology Facilities Council of the United Kingdom, whilst GI thanks the Royal Society for the award of a Newton Fellowship. The research leading to these results has received partial support from the European Research Council under the European Union’s Seventh Framework Programme (FP/2007-2013) / ERC Grant Agreement no. 308024. 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 the DPAC has been provided by national institutions, in particular the institutions participating in the Gaia Multilateral Agreement.