The Astrophysical Journal, 946:118 (17pp), 2023 April 1 https://doi.org/10.3847/1538-4357/acbb00
© 2023. The Author(s). Published by the American Astronomical Society.
The Imprint of Clump Formation at High Redshift. II. The Chemistry of the Bulge
Victor P. Debattista1 , David J. Liddicott1, Oscar A. Gonzalez2 , Leandro Beraldo e Silva1,3 , João A. S. Amarante4,1,16 ,
Ilin Lazar5, Manuela Zoccali6,7 , Elena Valenti8,9 , Deanne B. Fisher10 , Tigran Khachaturyants1 , David L. Nidever11 ,
Thomas R. Quinn12 , Min Du13 , and Susan Kassin14,15
1 Jeremiah Horrocks Institute, University of Central Lancashire, Preston, PR1 2HE, UK
2 UK Astronomy Technology Centre, Royal Observatory, Blackford Hill, Edinburgh, EH9 3HJ, UK
3 Department of Astronomy, University of Michigan, 1085 S. University Avenue, Ann Arbor, MI 48109, USA
4 Institut de Ciencies del Cosmos (ICCUB), Universitat de Barcelona (IEEC-UB), Martí i Franquès 1, E-08028 Barcelona, Spain
5 Centre for Astrophysics Research, School of Physics, Astronomy and Mathematics, University of Hertfordshire, Hatfield AL10 9AB, UK
6 Instituto de Astrofísica, Pontificia Universidad Católica de Chile, Av. Vicuña Mackenna 4860, 782-0436 Macul, Santiago, Chile
7 Millennium Institute of Astrophysics, Av. Vicuña Mackenna 4860, 82-0436 Macul, Santiago, Chile
8 European Southern Observatory, Karl Schwarzschild-Straße 2, D-85748 Garching bei München, Germany
9 Excellence Cluster ORIGINS, Boltzmann-Straße 2, D-85748 Garching bei München, Germany
10 Centre for Astrophysics and Supercomputing, Swinburne University of Technology, P.O. Box 218, Hawthorn, VIC 3122, Australia
11 Department of Physics, Montana State University, P.O. Box 173840, Bozeman, MT 59717, USA
12 Astronomy Department, University of Washington, Box 351580, Seattle, WA 98195, USA
13 Department of Astronomy, Xiamen University, Xiamen, Fujian 361005, Peopleʼs Republic of China
14 Space Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
15 Johns Hopkins University, 3400 North Charles Street, Baltimore, MD 21218, USA
Received 2022 August 29; revised 2023 January 11; accepted 2023 February 8; published 2023 April 7
Abstract
In Paper I, we showed that clumps in high-redshift galaxies, having a high star formation rate density (ΣSFR),
produce disks with two tracks in the [Fe/H]–[α/Fe] chemical space, similar to that of the Milky Way’s (MW’s)
thin+thick disks. Here we investigate the effect of clumps on the bulge’s chemistry. The chemistry of the MW’s
bulge is comprised of a single track with two density peaks separated by a trough. We show that the bulge
chemistry of an N-body + smoothed particle hydrodynamics clumpy simulation also has a single track. Star
formation within the bulge is itself in the high-ΣSFR clumpy mode, which ensures that the bulge’s chemical track
follows that of the thick disk at low [Fe/H] and then extends to high [Fe/H], where it peaks. The peak at low
metallicity instead is comprised of a mixture of in situ stars and stars accreted via clumps. As a result, the trough
between the peaks occurs at the end of the thick disk track. We find that the high-metallicity peak dominates near
the mid-plane and declines in relative importance with height, as in the MW. The bulge is already rapidly rotating
by the end of the clump epoch, with higher rotation at low [α/Fe]. Thus clumpy star formation is able to
simultaneously explain the chemodynamic trends of the MW’s bulge, thin+thick disks, and the splash.
Unified Astronomy Thesaurus concepts: Galactic bulge (2041); Milky Way formation (1053); Milky Way
evolution (1052); Milky Way dynamics (1051); Galaxy bulges (578)
1. Introduction with two peaks and a trough between them. In contrast, in the
Solar Neighborhood, two tracks17 are evident: at fixed [Fe/H],
The chemistry of the Milky Way’s (MW) bulge provides a high-[α/Fe] track corresponds to the thick disk, and a low-
important clues about its formation. The early measurements of [α/Fe] track corresponds to the thin disk. The bulge chemistry
Rich (1988), McWilliam & Rich (1994) established that the follows the thick disk track at low metallicity (Melendez et al.
bulge’s metallicity distribution function (MDF) is broad, 2008; Alves-Brito et al. 2010; Bensby et al. 2010; Hill et al.
reaching supersolar metallicities. More recent observations 2011; Bensby et al. 2013), but then extends to the most metal-
have shown that the MDF is at least bimodal with possible rich thin-disk stars. The location of the knee in the [Fe/H]–[α/
hints of additional peaks (Ness et al. 2013; Schultheis et al. Fe] plane has generally been found to be identical between the
2017; Rojas-Arriagada et al. 2020; Johnson et al. 2022), bulge and thick disk (Jonsson et al. 2017; Zasowski et al.
although this may partly be due to fitting multiple Gaussians to 2019), with perhaps minor differences (Johnson et al. 2014;
an intrinsically skewed distribution. Spectroscopic surveys Bensby et al. 2017; Schultheis et al. 2017), which may be
such as ARGOS (Freeman et al. 2013), GIBS (Zoccali et al. partly attributed to comparing bulge giants with local thick disk
2014), and APOGEE (Majewski et al. 2016) have mapped the dwarfs. Williams et al. (2016) found bimodalities in the bulge’s
chemistry across the bulge (e.g., Ness et al. 2013; Gonzalez [Fe/H] and [α/Fe] in the Gaia-ESO data, with the metal-rich
et al. 2015; Zoccali et al. 2017; Queiroz et al. 2021) generally stars exhibiting lower velocity dispersions than the metal-poor
finding that its [Fe/H]–[α/Fe] plane exhibits a single track, ones. The advent of the large APOGEE DR17 data set, and
matching data from Gaia Data Release 2 (DR2), has permitted
16 UCLan Visiting Fellow. more detailed studies of the bulge chemistry. Lian et al. (2020)
used the bulge’s chemistry to model its star formation history
Original content from this work may be used under the terms
of the Creative Commons Attribution 4.0 licence. Any further
distribution of this work must maintain attribution to the author(s) and the title 17 Different authors prefer either the term tracks or sequences to refer to the
of the work, journal citation and DOI. same thing. Throughout we will refer to tracks.
1
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
(SFH) and concluded that it is comprised of three phases: an presented in Section 3. We discuss our results, and give a brief
early high star formation rate (SFR) phase, which is interrupted summary of the main results, in Section 4.
by a quenched phase, which produces a gap in the chemistry,
followed by a later secular phase of low SFR. 2. The Simulations
The chemistry of the disk(s) differs from these trends. Many
explanations have been advanced for the disk α-bimodality. We use the clumpy simulation of Paper I, as well as a control
The “two-infall” model of Chiappini et al. (1997; see also simulation that fails to produce long-lived clumps; both these
Chiappini 2009; Bekki & Tsujimoto 2011; Tsujimoto & models are described in Beraldo e Silva et al. (2020). The
Bekki 2012; Grisoni et al. 2017; Khoperskov et al. 2021; subgrid physics of the two models differs only in the strength
Spitoni et al. 2021) suggests that a high SFR episode formed of the feedback employed. Both models are evolved from the
the high-α track, followed, around 8 Gyr ago, by a drop in the same initial conditions, comprised of a cospatial hot gas corona
SFR and then the infall of pristine gas that diluted the overall and dark matter halo with Navarro–Frenk–White (Navarro
metallicity of the MW, giving rise to the low-α population. et al. 1997) profiles. The dark matter halo has virial mass of
12
Recent work has focused on forming multiple chemical tracks 10 Me and a virial radius r200; 200 kpc. The gas corona,
via some variant of accretion events (Snaith et al. 2016; Grand which constitutes 10% of the mass within the virial radius,
et al. 2017; Mackereth et al. 2018; Buck 2020), including those starts with spin λ= 0.065 (Bullock et al. 2001), and as it cools,
of stars born out of the plane of the disk (Agertz et al. 2021). via metal line cooling (Shen et al. 2010b), it settles into a disk.
−3
In Clarke et al. (2019, hereafter Paper I), we presented a Stars form from dense gas (density >1cm ) when the
simulation of an isolated galaxy that produced a disk chemical temperature drops below 15,000 K, and the flow is convergent.
Gas particles are not allowed to cool below the resolution limit
dichotomy similar to the MW’s chemical thin+thick disks. At
( by setting a pressure floor p = 3Gò2ρ2, where G is Newton’s
early times largely over the first 2 Gyr, but continuing to 4 Gyr floor
gravitational constant, ò is the softening length, set at 50 pc, and
at a lower rate), the model develops clumps with high SFR
ρ is the gas particle’s density (Agertz et al. 2009). The feedback
densities, ΣSFR. The masses and SFRs of the clumps in this via supernovae Types Ia and II uses the blastwave prescription
model are comparable to those observed in high-redshift of Stinson et al. (2006). In the clumpy model, we couple 10%
galaxies (e.g., Guo et al. 2015; Dessauges-Zavadsky et al. of the 1051 erg per supernova to the interstellar medium as
2017; Cava et al. 2018; Guo et al. 2018; Huertas-Company thermal energy. In contrast, in the high-feedback model, 80%
et al. 2020). The clumps represent a second mode of star of the feedback energy is coupled to the gas. As shown in
formation, separate from the usual distributed star formation, previous studies (Genel et al. 2012; Hopkins et al. 2012; Buck
with high ΣSFR, leading to two tracks in the chemical, [Fe/H]– et al. 2017; Oklopčić et al. 2017), high feedback coupling
[α/Fe], plane. The rate of clump formation declines rapidly as inhibits the clumps, and Beraldo e Silva et al. (2020) show that
the gas fraction drops, thereby resembling the two-infall model. in that case the geometric properties of the disk(s) do not
In agreement with Bournaud et al. (2009), Paper I showed that resemble those of the MW. Feedback via asymptotic giant
clumps produce a geometric thick disk. The chemical and branch stars is also included. Gas chemical and thermal
geometric properties of the thick disk formed this way are diffusion uses the method of Shen et al. (2010b).
consistent with those of the MW (Beraldo e Silva et al. 2020). We evolve the models in isolation using a smooth particle
Moreover, Amarante et al. (2020) showed that the resulting low hydrodynamics +N-body tree-code based on GASOLINE
angular momentum tail of the old stars is consistent with the (Wadsley et al. 2004). The initial models are comprised of
“Splash” population in the MW (Di Matteo et al. 2019b; 106 particles in both the dark matter and gas components; both
Belokurov et al. 2020). models form ∼2× 106 stars. The clumpy model forms clumps
Paper I showed that some of the clumps sink to the center of during the first 2 Gyr, continuing at a lower rate to 4 Gyr, as
the galaxy, where they contribute to the formation of a bulge. shown in Paper I. The final disk galaxy has a rotational velocity
While definitively determining if clumps are long lived enough of 242 km s−1 at the Solar Neighborhood, making it compar-
to build bulges is challenging due to observational systematics able to the MW (see Figure 2 of Paper I). The high-feedback
(see, for instance, the discussion in Bournaud et al. 2014), model evolves without forming any significant long-lived
observations of the stellar populations (e.g., Guo et al. 2018; clumps. Henceforth we refer to the two models as the clumpy
Lenkic et al. 2021) and gradients of clump mass (Huertas- and high-feedback models.
Company et al. 2020; Ambachew et al. 2022) suggest that at Neither of these two models forms a bar. The formation of a
least some fraction of clumps likely do survive long enough to bar quenches star formation within most of the body of the bar
fall into the bulge. The chemistry of bulges formed with a (e.g., Khoperskov et al. 2018). In order to compare with the
significant contribution from clumps has not been studied MW, we assume that the MW’s bar formed at t= 6 Gyr (which
extensively in the literature, despite frequent suggestions that would make it ∼8 Gyr old now).
bulges, including the MW’s, may be partly built from clumps
(e.g., Nataf 2017; Queiroz et al. 2021). Interestingly, Immeli 3. Bulge Stellar Populations
et al. (2004) found a bimodal distribution of [Mg/Fe] within
the bulge of their clumpy chemodynamical model. Inoue & 3.1. The Chemistry of the Bulge
Saitoh (2012) found a metal-rich bulge formed from clumps but The top left panel of Figure 1 presents the chemistry of the
did not study the chemistry in greater detail. Therefore in this stars within a galactocentric radius R= 1 kpc at 10 Gyr in the
paper we study the consequences of star formation in a clumpy clumpy model. As in Paper I, we apply Gaussian measurement
mode on the chemistry of the bulge. uncertainties of σ[Fe/H] = 0.1, and σ[O/Fe] = 0.03 to mimic
The paper is organized as follows. Section 2 presents the the measurement errors in APOGEE (Nidever et al. 2014). The
simulations used in this paper. The chemistry, star formation, chemical space has a single track, with the density peaked at
kinematics, and spatial variation of the model bulges are two locations: one metal-rich at [Fe/H]; 0.55 and a broader
2
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 1. The density of stars in the [Fe/H]–[O/Fe] chemical space at t = 10 Gyr. Top: all stars within R = 1 kpc. Bottom: stars at R > 5 kpc with 1000 random
bulge stars superposed as red points. At left is the clumpy model, while at right is the high-feedback model. Smoothing in [Fe/H] and [O/Fe] has been applied to all
panels to match the chemical resolution of APOGEE, as described in Section 3.1. Bins with less than 100 stars have been suppressed. In the clumpy model, the bulge
chemistry matches that of the disk in the high-[O/Fe] region, in agreement with MW trends, and in contrast to the high-feedback model.
metal-poor peak at [Fe/H];−0.1. The bottom left panel of The main difference is at early times, when the presence of the
Figure 1 presents the chemistry of the clumpy model’s thin clumps briefly raises the overall peak SFR by ∼20%. In the
+thick disks at R> 5 kpc, and compares this with the high-feedback model, these clumps are short-lived (Genel et al.
chemistry of the model’s bulge (the red points represent a 2012; Hopkins et al. 2012; Buck et al. 2017; Oklopčić et al.
random selection of 1000 bulge particles). The chemistry of the 2017), and the SFR is therefore briefly lower.
bulge follows that of the thick disk at [Fe/H] 0, and then
continues to more metal-rich than the thin disk. The MW’s
bulge exhibits the same trend (e.g., Melendez et al. 2008; 3.2. Evolution of the Bulge’s Chemistry
Alves-Brito et al. 2010; Bensby et al. 2010; Hill et al. 2011; The fact that the clumpy model’s bulge chemistry has a
Bensby et al. 2013; Lian et al. 2020). We have verified that the single, double-peaked track that matches that of the thick disk
trends in Figure 1 are already in place by t= 6 Gyr. at [Fe/H] 0 is strikingly similar to what is observed in the
The right panels of Figure 1 present the chemistry of the MW. Understanding this trend therefore can help unravel the
high-feedback model. A number of important differences formation of the MW’s bulge. Thus we next explore the
between the clumpy and high-feedback models are evident. evolution of the bulge chemistry to understand how clumpy
The first difference is that the track of the bulge in chemical star formation produces these trends.
space no longer has two peaks. Instead the bulge has a single Figure 3 shows the chemical evolution of the bulge inside
sharp peak at [Fe/H]; 0.6 with a long tail to lower R= 1 kpc for both models. We show the MDF and the α
metallicities. Moreover, this model does not have a bimodal distribution function (α DF) for all bulge stars formed up to 2,
chemical distribution in the disk (see also Beraldo e Silva et al. 4, 6, 8, and 10 Gyr. The clumpy model, at t= 4 Gyr, when
2020), which happens because the high-α stars form only via clump formation fully ceases, has a bulge MDF, which is
the clumpy star formation mode in these simulations. As a bimodal (top left panel), with a low-metallicity peak at [Fe/
consequence, the bulge chemical distribution is offset vertically H];−0.1 and a small peak at [Fe/H]; 0.4. The high-
in [O/Fe] relative to the disk. While the bulge has a high SFR metallicity peak grows in importance as subsequent in situ star
and can therefore reach a high [O/Fe], this is not the case in the formation adds a population of high-metallicity stars. The
disk, and the bulge ends up more α-rich than the disk. The lack trough between the two peaks falls at [Fe/H]; 0.25. In the
of a trough in the bulge’s chemistry and the difference between MW’s bulge, the metallicity of the trough varies with position
the bulge’s peak α and that of the disk are different from the in the range [Fe/H]∼−0.2 to 0.2 (Zoccali et al. 2017). After
trends observed in the MW. t= 4 Gyr, the α DF of the clumpy model (bottom left panel)
In spite of these differences in chemical space, the overall has a fixed peak at high [O/Fe] (at ≈0, but we caution that [O/
SFH of these two models is very similar, as seen in Figure 2. Fe] values often have significant offsets in simulations
3
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 2. The overall star formation history of the clumpy and high-feedback models.
compared to observations, as we also found in Paper I). At infall of clumps into the bulge of the clumpy model between 2
t= 4 Gyr, the α DF has a point of inflection at low [O/Fe], and 4 Gyr. After 4 Gyr, when no further clumps form in the
where a pronounced second peak later develops. A double- disk of the clumpy model, the chemical evolution of the two
peaked αDF is similarly present in the MW’s bulge (e.g., Lian bulges proceeds very similarly, with an increasing number of
et al. 2020). stars at the high-[Fe/H], low-[α/Fe] peaks. During this time,
In contrast, the chemical evolution of the high-feedback the clumpy model develops a second peak at low [α/Fe],
model (right panels) results in only a single peak in the bulge’s which had formed earlier in the high-feedback model. In
MDF, and only a weak double peak in the bulge α DF. At best summary, it is not the differences in their SFRs that give rise to
a weak trough is visible in chemistry of its bulge. The two the different chemistries of the two bulges, but the infall of
models differ at the low-[Fe/H] peak (i.e., at the high-[O/Fe] clumps onto the bulge of the clumpy model, which drives the
peak), which must represent the location where the clump continued growth of the low-[Fe/H] peak in its chemistry.
formation plays an important role in one model and is absent
from the other. 3.3. Formation Location
Small differences between the clumpy and the high-feedback
models are already present at 2 Gyr, which Figure 2 shows has Paper I showed that the clumps in the clumpy simulation
the largest differences between the global SFRs of the two often fall to the center. If clumps are disrupted before they can
models. At 2 Gyr, the MDF of the clumpy bulge has a peak at reach the bulge, then they may play a less prominent role in the
low [Fe/H], while a peak at high [Fe/H] is incipient, but not formation of the bulge. We therefore consider the formation
yet prominent. The high-feedback bulge has a very similar location of bulge stars to test the effect of the infalling clumps
MDF, but it has only a single peak at roughly the same subsolar on the chemistry of the bulge.
[Fe/H] as in the clumpy model. The low-[Fe/H] peak is more The top left panel of Figure 4 shows the distribution of the
prominent in the clumpy bulge than that in the high-feedback formation radius, áRformñ, in the chemical space. An important
bulge, but the overall trends are similar. Similarly the DFs of conclusion from this plot is the different origins of the two
α
the two models are not yet very different, with a single peak at MDF peaks. The stars at the low-[Fe/H] peak in the chemical
high α. The differences between the chemistry of the two track have large áRformñ, indicating that many of them are
bulges become larger between 2 and 4 Gyr, despite the fact that forming outside the bulge and reaching it via clumps. The high-
[Fe/H] peak instead is produced by in situ18 star formation (as
the global SFRs of the two models are more similar at these in the high-feedback model; seen in the top right panel of
times. In the clumpy model, the separate peak at high [Fe/H] Figure 3). The bottom left panel of Figure 4 shows the fraction
now becomes more developed, while the continuing enrich- of bulge stars that are ex situ. This is high at intermediate [Fe/
ment in the bulge of the high-feedback model results in only a H] and low at high [Fe/H], closely mirroring the top left panel.
single peak at high [Fe/H]. The low-[Fe/H] peak in the clumpy When we plot the distribution of stars in the bulge’s
model grows in importance at this time, while shifting to higher chemical space excluding those stars that formed outside
[Fe/H]. At the same metallicities as the low-[Fe/H] peak of the Rform= 2 kpc, which we do in the top right panel of Figure 4,
clumpy model, the bulge of the high-feedback model barely we find that the low-[Fe/H] peak is substantially reduced, with
changes during this time. In the high-feedback bulge, the α DF the track resembling somewhat the distribution of the high-
begins to develop a peak at low [α/Fe], while in the clumpy feedback model in Figure 1. (The small peaks remaining after
bulge the low-[α/Fe] peak has not yet started to be visible, but
the high-α peak continues to grow while shifting to lower [α/ 18 Here we use the terms in situ and ex situ to refer to formation inside or
Fe]. As we show below, the driver of these differences is the outside the bulge, but within the galaxy.
4
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 3. The evolution of the MDF (top) and α DF (bottom) of stars within the bulge (R „ 1 kpc) between t = 2 Gyr, and t = 10 Gyr. At left is the clumpy model,
and at right is the high-feedback one. In the clumpy model, a bimodality is present in the MDF at t = 10 Gyr with a broad, low peak at [Fe/H] ∼ −0.1 and a narrow,
high peak at [Fe/H] ∼ 0.5. The bimodality is already evident, although weaker, at t = 2 Gyr, when clump formation has started to die down, and is well established at
4 Gyr. A bimodality is also present in the α DF at t = 10 Gyr, with a broad, low peak at [O/Fe] ∼ 0 and a narrow, high peak at [O/Fe] ∼ −0.3. This bimodality is
significantly weaker and/or absent at t = 4 Gyr. In the high-feedback model, instead, only a single peak develops in the MDF although the α DF still has a weak
second peak. All distributions have been normalized to the corresponding peak at 10 Gyr. In the top row, the vertical dotted lines indicate the regions around the peaks
where we define MDF peaks discussed in Section 3.5.
this subtraction are caused by stars formed in clumps that are 3.4. The Link between the Single Track and the Star
still star-forming inside R= 2 kpc.) This explains why the Formation Mode
bulge chemistry at low [Fe/H] is such a good match to the In the left panels of Figure 5, we plot the density of stars in
chemistry of the thick disk: many of these stars share a similar the space of final versus formation radii (Rfinal versus Rform).
origin. The stars that form during the clump epoch, tform< 4 Gyr, and
The bottom right panel of Figure 4 shows the distribution of that end within the inner 1 kpc (top left panel) form at a range
those stars excluded from the second panel, i.e., the stars that of radii, including a significant contribution forming in situ. For
end within the bulge that formed at Rform> 2 kpc. This shows the stars formed after the clump epoch, 4„ tform/Gyr„ 6,
that the bulk of these ex situ stars arriving within clumps settle (bottom left panel) star formation occurs in situ, resulting in the
along the bulge track, with their highest density at the location diagonal distribution in the Rfinal–Rform space. In the absence of
of the low-metallicity peak. clumps, stars only reach the bulge from farther out via eccentric
A number of additional conclusions can be drawn from the orbits; the bottom left panel shows that the fraction of such
top left panel of Figure 4. First is the fact that clumps bring stars is low.
with them a small population of low-α stars, which settle below The stars in the bulge therefore are a mix of those formed
the low-[Fe/H] peak around [Fe/H]∼−0.4 and [O/ in situ and those accreted in clumps. Paper I showed that there
Fe]∼−0.1. As shown in Paper I (in Figures 15 and 17), some are two modes of star formation: a high ΣSFR and a low ΣSFR
of the stars formed in clumps have low [α/Fe]. This population one. Clumps are associated with high ΣSFR (ΣSFR
of stars is relatively small and does not contaminate the 1 Me yr−1 kpc−2) while lower ΣSFR is typical of distributed
chemical distribution significantly. The second point is that, to (nonclumpy) star formation (see Figure 15 of Paper I). The
a large extent, the chemistry of the bulge, at the high- and low- right panels of Figure 5 show the distribution of áSSFRñ in the
metallicity ends, is dominated by stars formed in situ, and is same Rfinal versus Rform space. For stars with tform< 4 Gyr (top
contaminated by clumps only at −0.5 [Fe/H] 0.0. right panel), the high áSSFRñ at Rfinal< 0.5 kpc is produced by
5
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 4. Top left: distribution of áRformñ in the chemical space of stars formed in the first 4 Gyr that end within the inner 1 kpc of the clumpy model. The accreted
clumps are responsible for the low-[Fe/H] peak while in situ star formation produces the high-[Fe/H] peak. The contours indicate the density of particles; the 5
contour levels span a factor of 10. Bottom left: the fraction of ex situ stars (those with Rform > 2 kpc) that end up in the bulge (Rfinal < 1 kpc). Top right: the in situ
bulge, showing the distribution of stars contained within Rfinal < 1 kpc when stars with Rform > 2 kpc are excluded. While the distribution is not completely smooth,
no prominent peak at low-[Fe/H] is evident. Bottom right: the ex situ bulge, defined as those stars within Rfinal < 1 kpc with Rform > 2 kpc.
the full range of Rform, which therefore must include stars Why then does the in situ star formation not produce a
formed in clumps that have fallen in, as well as those formed separate track in the bulge chemistry like the disk’s low-[α/Fe]
in situ. This is made clearer by comparing with stars that form track? The bottom panel of Figure 6 plots áSSFRñ for stars that
after 4 Gyr (bottom right panel) when clump formation has were born within a given radius by tform= 4 Gyr. The vast
ceased; now bulge stars have relatively low áSSFRñ  majority of early stars formed within 2 kpc formed via the
3M yr-1 kpc-2 (except at the very center), and high-ΣSFR mode. Therefore, the early bulge itself acts as a

R clump of high Σ , as first shown by Mandelker et al. (2014).
form; Rfinal. These stars clearly are forming in situ rather SFR
than falling in as clumps. Thus the bulge never gets to form a low-α track: even in the
Therefore stars born outside the bulge in the clumpy absence of clumps falling into the bulge, for instance because
high-Σ mode are reaching the bulge. In the top panel of they are disrupted before they reach the center, the bulge
SFR
Figure 6, we plot Σ for stars forming before 4 Gyr that end chemistry will still lack metal-poor low-α stars. Indeed even
SFR
at different radii within the inner galaxy. From R „ 2 kpc the high-feedback model has only a single track in the bulge,
final
(red curve) to Rfinal„ 1 kpc (green curve), the contribution of and is α-rich (compared with the disk), as can be seen in the top
right panel of Figure 1.
the high-ΣSFR mode of star formation rises, and overwhel-
mingly dominates at Rfinal< 0.5 kpc (blue curve). This would
seem to imply that infalling clumps dominate the inner galaxy. 3.5. Bulge Ages and Quenching
However, the middle panel of Figure 6, which shows the If star formation continues in the bulge after the clump epoch
cumulative distribution of Rform for stars at different Rfinal, ends, then the younger stars will necessarily be at the high-[Fe/
shows that a significant in situ population is also present. H] peak. The high-[Fe/H] peak then would be younger, on
Indeed the fraction of stars that formed within Rform= 2 kpc average, than the low-[Fe/H] peak. In the MW, the difference
rises as Rfinal decreases, although it never exceeds ∼60%. in mean age between the two peaks would be governed by
Roughly half the stars that end up at Rfinal„ 2 kpc were born when the bar forms, because bars generally quench star
outside this region. Thus clumps are delivering a significant formation within most of their radius (including the vertically
fraction of the bulge’s mass, but in situ star formation is equally thickened part that forms the bulge). Figure 7 shows the SFH
important. up to 6 Gyr for the stars that end within R= 1 kpc (we have
6
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 5. The number of stars (left column) and áSSFRñ in the Rform–Rfinal space for stars at the center of the clumpy model. The top row shows the distributions for
tform „ 4 Gyr while the bottom row is for 4 < tform/ Gyr „ 6. The diagonal structure in the bottom panels shows predominantly in situ star formation after
tform = 4 Gyr, whereas the upper panels show a significant population of bulge stars brought in by clumps. (Note the different scales on the two axes. The diagonal
dashed green line indicates Rform = Rfinal.)
checked that the result does not change qualitatively if we chemical track (as can be seen in Figure 1), suggesting that the
consider stars inside R= 2 kpc). This shows that the high- SFH need not be responsible for the trough.
metallicity population ([Fe/H]… 0.25) overlaps in age with the The evolution of the bulge MDF in the clumpy model, seen
low-[Fe/H] population. However no new stars with low [Fe/ in Figure 3, shows that the bulge reaches the high-[Fe/H]
H] form after ∼4 Gyr. regime already by 2 Gyr, and is bimodal already by that point.
Lian et al. (2020) interpreted the trough between the two The bimodality increases at later times, particularly after 4 Gyr,
peaks in the MW bulge’s chemical track as an episode of but the trough is present before the clumpy episode is over.
quenching in its SFH, before star formation restarted and Thus the high-[Fe/H] peak contains old stars and represents the
produced the high-[Fe/H] peak. The top panel of Figure 7 ordinary chemical evolution of a rapidly star-forming system. If
shows that the star formation in the clumpy model’s bulge we understand the [Fe/H]-enrichment as developing smoothly,
never drops to zero, although its chemical track in Figure 1 then stars in the trough will be, on average, slightly older than
develops a trough between the two peaks. The bottom panel of those in the high-[Fe/H] peak. The stars at the low-[Fe/H]
Figure 7 presents the SFH of the high-feedback model. Despite peak, because they are a mix of in situ stars and stars accreted
the similarity in the SFH of the two models, only the clumpy via clumps, represent a range of ages, from older than the
model develops a trough between two peaks in the bulge trough (from the in situ evolution) to stars younger than the old
7
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 7. Star formation history to 6 Gyr, i.e., 2 Gyr after the end of the clump
era in the clumpy model, for stars that end up within Rfinal < 1 kpc, separated
by metallicity. The full distribution is shown in black. The separation into
metal-rich and metal-poor is at [Fe/H] = 0.25, which marks the minimum
between the two peaks in the MDF of the clumpy model. The stars in the low-
[Fe/H] and high-[Fe/H] populations are shown in blue and red, respectively.
The top panel shows the clumpy model while the bottom panel shows the high-
feedback model split at the same metallicity. The overall similarity of the SFHs
of the two models suggests that quenching is not responsible for the trough in
the chemical space of the clumpy model, seen in the top left panel of Figure 1.
region where the low-[Fe/H] peak would be has predominantly
older stars and only overlaps the high-[Fe/H] peak’s ages in
the exponential wing of the distribution.
We conclude that the trough between the high- and low-
Figure 6. Star formation modes at the center of the clumpy model for stars with
t „ 4 Gyr. Top: the normalized distribution of the star formation rate [Fe/H] peaks in the clumpy model is not due to a quenching
form
density for different final radii. Middle: the cumulative formation radius of stars of in situ star formation but rather due to the end of clumps
that end up at different radii. Bottom: the normalized distribution of star delivering stars to the low-[Fe/H] peak of the bulge. A
formation rate density for different formation radii. In the top and bottom possible diagnostic of this scenario is the age distribution of
panels, the black histograms refer to all the stars within the model. All
distributions use kernel density estimates with a Gaussian kernel and window the low-[Fe/H] peak compared with that in the trough: stars at
width satisfying Silverman’s rule (Silverman 1986). the low-[Fe/H] peak should include younger stars than those
in the trough. We explore this prediction for the clumpy
stars in the high-[Fe/H] peak (from the later stages of clump model in the top panel of Figure 9, where we plot the mean
time of formation, átformñ, of stars in the chemical space of the
accretion). Figure 8 shows the distribution of ages at the two bulge stars formed by t= 6 Gyr. Along the ridge of the
peaks, within the [Fe/H] limits indicated by the vertical dotted
 chemical track from metal-poor to metal-rich, we reach a local
lines in the MDFs of Figure 3. At 2 tform/Gyr 4, the ages maximum in átformñ at the location of the low-[Fe/H] peak,
of stars at the high-[Fe/H] peak significantly overlap those of while the high-[Fe/H] peak is the location of late star
the youngest stars at the low-[Fe/H] peak. We show, as dashed formation and has the largest átformñ (i.e., youngest stars). In
lines, the age distributions for the high-feedback model in the between, at the trough, átformñ has a local minimum, meaning
same metallicity ranges. While the age distribution at the high- the stars in this region are older. Observationally, this dip
[Fe/H] peak is comparable to that of the clumpy model, the gives the appearance of a drop in the SFR of the bulge and
8
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 8. The SFH of stars centered at the clumpy model’s two MDF peaks within the inner 1 kpc. The clumpy model is shown by the solid lines while the high-
feedback model at the same [Fe/H] ranges is shown by the dashed lines. The two metallicity ranges selected are indicated by the vertical dotted lines in the top row of
Figure 3.
thus resembles a quenching episode. However this is clearly anticipate that stars at the high-[Fe/H] peak would dominate near
not the case in the evolution of the clumpy model. In contrast, the mid-plane while the large heights would be dominated by the
the bottom panel shows the mean age of the high-feedback low-[Fe/H] peak if a bar had formed.
model, which shows that the mean age increases mono-
tonically along the ridge in this case.
A final noteworthy property of the clumpy bulge’s chemistry 3.7. Comparison with the Milky Way
is that the trough occurs just beyond the highest metallicity of Despite the absence of a bar in the clumpy model, we can
the thick disk track (see Figure 1). This happens because the compare the vertical distribution of the MDF with the MW’s.
trough is not significantly polluted by stars formed in the same In order to do this, we use the model at 6 Gyr assuming that the
clumps that produced the thick disk. bulge is quenched by bar formation at this time. We apply a
coordinate transformation of the model’s Cartesian coordinates
to Galactic coordinates, (l, b, d), after placing the Galactic
3.6. Dependence of Kinematics on Chemistry
center at 8 kpc from the Sun. We select particles across
The top panel of Figure 10 shows profiles of the radial constant longitude stripes (−6°.5< l< 6°.5) at different lati-
velocity dispersion, σR, of stars at t= 4 Gyr separated into [O/ tudes (1°.5< |b|< 2°.5 and 5°.5< |b|< 6°.5, restricted to a
Fe] bins. The high-[O/Fe] stars are generally hotter, by distance 7< d/ kpc< 9. This represents a typical spectroscopic
20–30 km s−1, than the low-[O/Fe] stars even just at the end of selection of giant stars in the MW bulge (e.g., Wylie et al.
the clumpy epoch. This reflects on the chaotic interaction of 2021) with which variations as a function of latitude are
clumps near the center of the galaxy, which heats the high-[O/ studied. Figure 12 shows the resulting distribution of the
Fe] populations at birth. selected stars in chemical space; these display two over-
The bottom panel of Figure 10 shows profiles of the mean densities that change their relative contributions as a function
streaming velocity, áVfñ; although clumps are falling to the of Galactic latitude, as in the MW.
center, the bulge remains rotationally supported, because the The chemical track in Figure 12 is comprised of a sequence
clumps are on in-plane, prograde orbits, which are known to of [Fe/H]-poor stars, followed by a trough and then a shorter
produce rapidly rotating remnants even when the resulting sequence of [Fe/H]-rich stars whose relative contribution
mergers are collisionless (e.g., Read et al. 2008; Hartmann decreases with increasing height from the plane. Without any
et al. 2011). In cosmological simulations, Inoue & Saitoh scaling applied to the simulation, the [Fe/H]-rich population is
(2012) also found rapidly rotating bulges forming from clumps. no longer present at a latitude of b= 6°. Since the simulated
The low-[O/Fe] stars are more rapidly rotating, by galaxy has not formed a bar, the detailed properties of the two
∼50–100 km s−1, as expected given their lower velocity populations and their spatial variations are not directly
dispersion. comparable to those in the MW. The specific distributions
We measure the actions of stars using AGAMA (Vasiliev 2019), seen in the MW would depend on many details, such as the
which uses the Stäckel fudge of Binney (2012), assuming a epoch of bar formation, the vertical thickening of the bar, and
flattened axisymmetric potential for the disk and a spherical the star and clump formation histories. However, the presence
potential for the halo. Figure 11 shows the mean radial action, of this overall bimodality and trend in the simulation is
áJRñ, in the chemical space, for stars in the inner 1 kpc at 6 Gyr. consistent with the observations of Rojas-Arriagada et al.
Bearing in mind that Debattista et al. (2020) found that bar (2019) and Wylie et al. (2021) based on [Mg/Fe] abundances
formation substantially steepens the vertical gradient of áJRñ, we from APOGEE and ARGOS data. They showed, from a large
9
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 9. The distribution of átformñ in the chemical space of stars formed in the first 6 Gyr that end within the inner 1 kpc of the clumpy (top) and high-feedback
(bottom) models. The location of the trough in the chemical space of the clumpy model corresponds to a local minimum in átformñ. The contours indicate the density of
particles; the 5 contour levels span a factor of 10.
number of stars, that the [Fe/H] bimodality is produced by a Despite not having formed a bar, the simulation trends
low-α sequence of stars over a range of ∼0.5 dex around a suggest that the chemical bimodality in the bulge produced by
solar metallicity that merges with the main high-α sequence. In clumps is plausibly able to account for the spatial variations of
the simulation, there is a third, much smaller component in the chemistry of the MW’s bulge.
Figure 12 that appears as a lower-[O/Fe] overdensity in the
metal-poor regime. This population becomes more important at 4. Discussion
higher distances from the plane but clearly remains a minor
component throughout. This component is comprised of stars The single track in the chemical space of the bulge, i.e., the
formed in clumps that have lower [O/Fe] that form at large absence of an [α/Fe] bimodality for a fixed [Fe/H], in the
radii (see Figure 4). Their appearance suggests that the clumpy clumpy simulation is similar to that observed in the MW,
including the fact that it has two density peaks along the track.
simulation underestimates the evaporation rate of the clumps, Together with the clumpy model’s two tracks in the chemical
possibly because the feedback should be higher. We note, space of the disk (Paper I), this is a striking agreement with the
however, that a hint of such a low-[O/Fe], metal-poor trends seen in the MW, and suggests that the model is capturing
population can be seen in Wylie et al. (2021) where an a generic behavior. Altogether, these results demonstrate that a
increasing width of the low-[Fe/H] sequence as a function of holistic view of the chemistry of the entire MW (both bulge and
height is evident in their Figure 25, but further studies with thin+thick disks) provides a more stringent constraint on how
higher number statistics are needed to confirm this. the early MW formed (see also Di Matteo 2016).
10
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Mackereth et al. 2018; Buck 2020). Mackereth et al. (2018)
found that such outcomes only occur in about 5% of the
EAGLE simulation galaxies, but Buck (2020) found them to be
more common in the NIHAO simulation suite. Likewise, the
Splash has been interpreted as a combination of accreted
material and the kinematically heated disk after the GSE
merger (Di Matteo et al. 2019a; Gallart et al. 2019; Mackereth
et al. 2019; Belokurov et al. 2020).
The bulge is an important test of the hypothesis that the GSE
merger is exclusively responsible for the chemodynamics of the
MW because, on the one hand, the merger cannot leave a
classical bulge more massive than ∼8% of the total stellar mass
(Shen et al. 2010a; Bland-Hawthorn & Gerhard 2016;
Debattista et al. 2017), while at the same time producing a
bulge chemistry with a single track with two peaks. To date,
cosmological simulations that produce the chemical thin and
thick disks appear to produce two, or more, tracks in the
chemistry of the bulge (e.g., Grand et al. 2018; Buck 2020). We
have shown here that the chemistry of the bulge can largely,
and very naturally, be produced by clumpy star formation
(including within the bulge itself). Thus most of the
chemodynamics of the early MW, excluding the accreted halo,
can now be explained by clumps. However, because the GSE
merger certainly happened, it is important to understand to
what extent a merger in the presence of clumps is able to
explain the details of the MW’s chemodynamics. We will be
exploring exactly this with project GASTRO (Amarante et al.
2022).
Further complicating matters, besides the GSE, there have
been suggestions of at least one other equally massive merger in
Figure 10. Radial velocity dispersions, , (top) and mean rotational velocity, the MW during its early evolution (Massari et al. 2019;
σR
áVfñ, (bottom) of stars at 4 Gyr, in the inner 2 kpc of the clumpy model, as a Forbes 2020; Horta et al. 2021). Horta et al. (2021) used
function of [O/Fe]. Low-[O/Fe] stars have lower radial velocity dispersions APOGEE DR16 and Gaia DR2 to characterize the stars of this
and higher rotation. merger event, which they called Heracles. 20
“ ” They estimated its
stellar mass as ∼5× 108 Me, i.e., as massive as GSE (see also
4.1. Comparison with Other Scenarios Kruijssen et al. 2020). The stars associated with Heracles are
We have shown that an episode of star formation in clumps located at R< 5 kpc, and are thus more bound to the Galactic
is able to explain the twin peaks in the bulge s single track in potential than the GSE remnant (but, see also Lane et al. 2022,
’
for a discussion of whether Heracles could be an artifact in the
the chemical space. The bulge track follows that of the thick
E–L plane of APOGEE’s selection function). These stars are
disk at low metallicity but then continues to the thin disk and z
also chemically distinct from the GSE (Horta et al. 2021; Naidu
beyond at high metallicity, as observed in the MW. In Paper I et al. 2022) and would imprint as bursts in the SFH of the inner
and Beraldo e Silva et al. (2020), we showed that the chemical MW (Orkney et al. 2022). Naidu et al. (2022) estimated it was
thick and thin disks produced via clumps have similar accreted ∼1.7 Gyr before GSE. Recently, Myeong et al. (2022)
properties (scale-lengths and scale-heights, kinematics, and argued for an in situ origin of Heracles. More recently, this
MDFs) as found in the MW. Amarante et al. (2020) showed population has been interpreted as the first stars that formed in
that clumps also produce the relatively metal-rich population the MW, based on Gaia, APOGEE DR17, and H3 survey data
that bridges the thick disk and the inner halo, which has been (Belokurov & Kravtsov 2022; Conroy et al. 2022; Rix et al.
termed “the Splash” (Di Matteo et al. 2019a; Belokurov et al.
19 2022).
2020). The clump scenario predicts that the thin and thick An alternative, popular model for the formation of the
disks were forming at the same time (Paper I), which appears to geometric thick disk posits that it formed in situ, already thick,
be consistent with the presence of RR Lyrae in the thin disk as in an “upside-down, inside-out” manner (Bird et al.
well as an age overlap between the chemical thin and thick 2013, 2021). Support for this model includes the short scale-
disks (Beraldo e Silva et al. 2021). length of the (chemical) thick disk (Bovy et al. 2012; Hayden
Since Gaia’s confirmation of the Gaia-Sausage-Enceladus et al. 2015) and the homogeneity of the high-α population. This
(hereafter GSE; Belokurov et al. 2018; Helmi et al. 2018) scenario is also supported by the high gas velocity dispersions
merger remnant, the chemodynamics of the early MW have and star formation rates in high-redshift galaxies (e.g., Kassin
been interpreted as products solely of this merger. Numerous et al. 2012; Wisnioski et al. 2015). The lack of flaring in the
cosmological simulations have indeed shown that disk high-α populations is also consistent with the upside-down
chemical bimodalities can arise from gas-rich mergers (e.g., scenario (Bovy et al. 2016; Ted Mackereth et al. 2017; but see
Brook et al. 2005; Snaith et al. 2016; Grand et al. 2018;
20 Massari et al. (2019), Forbes (2020) dubbed this remnant “Kraken” and
19 Di Matteo et al. (2019a) refer to this feature as “The Plume.” “Koala,” respectively.
11
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 11. The mean radial action, áJRñ, in the chemical space of the clumpy model. The contours indicate the density of particles; the 5 contour levels span a factor
of 10.
also Lian et al. 2022). In general, however, studies of the during the first infall, producing the high-α sequence, followed
upside-down formation scenario have provided no explanation by a second infall with low SFR, producing the low-α
for the disk chemical bimodality, or the chemistry of the bulge. sequence. As noted in Paper I, the clump model is similar, in
An alternative flavor of the upside-down formation scenario is terms of SFR, to this model, and the outcome may be
based on misaligned star formation. Meng & Gnedin (2021) indistinguishable. However, our results for the nonclumpy
showed that, in their cosmological simulations, stars always model, which has an SFH not much different from that of the
form in a thin disk, even at z> 1.5, and only give the clumpy model, but which fails to form a disk chemical
appearance of an upside-down formation because disks tilt bimodality, are at odds with a pure early high SFR producing a
rapidly at early times, which leaves the star-forming plane disk chemical bimodality. Clumps produce the chemical
misaligned (warped) with respect to the main disk plane. The bimodality by boosting the star formation rate density by a
subsequent precession of the stars formed off the plane factor of ∼100 compared to distributed star formation (Clarke
continuously inflates the height of the main disk (see also et al. 2019). Khoperskov et al. (2021) presented several
Khachaturyants et al. 2021). More recently, Tamfal et al. simulations that produced a thin+thick disk chemical bimod-
(2022) used a high-resolution (∼109 particles) zoom-in ality, which they attributed to the rapidly dropping SFR,
cosmological simulation to show that the disk is already coupled with outflows (see also Vincenzo & Kobayashi 2020)
forming thin as early as z∼ 7–8, with no upside-down similar to the two-infall model. The authors also noted that their
formation. This rotationally supported disk thickens slowly models undergo a period of clump formation, with comparable
due to internal instabilities and external perturbations, with clump masses to what we found in Paper I.
stellar accretion from satellites providing the main geometric
thick disk. In a similar vein, Agertz et al. (2021; see also 4.2. Observational Tests
Renaud et al. 2021a, 2021b) proposed that the origin of the Clumps are observed in more than half of high-redshift MW
chemical bimodality of the thin+thick disks is due to different progenitors (e.g., Elmegreen & Elmegreen 2005; Ravindranath
chemistry in the inner and outer disks, which accreted their gas et al. 2006; Elmegreen et al. 2007; Forster Schreiber et al.
from separate filaments. Early rapid star formation and mergers 2011; Genzel et al. 2011; Guo et al. 2012, 2015). Observed at
in the inner disk gave rise to the high-α thick disk population, high resolution, clumps are found to have sizes of order
while star formation in the outer misaligned disk is inhibited by 100–500 pc, average masses of ∼108Me (Livermore et al.
the low density of the gas until the last major merger triggers 2012, 2015; Fisher et al. 2017; Cava et al. 2018) and contribute
star formation in the outer disk, which becomes the metal-poor, about 7% of the ongoing star formation rate (Wuyts et al.
low-α thin disk. The continuing star formation then builds the 2012). Aside from the formation of a geometric thick disk
metal-rich, low-α thin disk we see today. Renaud et al. (2021b) (Bournaud et al. 2009; Clarke et al. 2019; Beraldo e Silva et al.
presented the chemistry of this simulation; the bimodal tracks 2020), the presence of clumps does not lead to substantial
extend to the inner galaxy, contrary to what is observed in the differences in the morphological properties of galaxies. Indeed
MW. It is unclear whether this outcome can be avoided in this the cosmological zoom-in simulations of Inoue & Yoshida
scenario. (2019), with identical initial conditions but varying gas physics,
The classical two-infall scenario of Chiappini et al. (1997; found a strong dependence of clump formation on the equation
see also Chiappini 2009; Bekki & Tsujimoto 2011; Tsujimoto of state of the gas, but very little effect on the global properties
& Bekki 2012; Grisoni et al. 2017; Spitoni et al. 2021) of the galaxies. The signatures of clumps are therefore
proposes that the formation of two sequences in the disk primarily chemical, because the masses of the clumps
chemistry results from two-infall episodes, with high SFR are modest (Livermore et al. 2012; Fisher et al. 2017;
12
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Figure 12. Gaussian kernel estimate of the density in the [Fe/H]–[O/Fe] chemical space for stars in the clumpy model. Ten equally spaced contours show the density
distribution of stars selected to be within the volume bounded by −6°. 5 < l < 6°. 5 and 7 < d/ kpc < 9, in 1° stripes at b = 2° (bottom) and 6° (top). The histogram at
the top shows the full distribution of [Fe/H]; the minimum between the two [Fe/H] peaks, indicated by the dotted line in the central panel, splits the distribution into
the low- and high-[Fe/H] populations. The [O/Fe] histograms of these two populations are shown at the right with the full (black), low [Fe/H] (red), and high [Fe/H]
(dotted blue).
Cava et al. 2018; Benincasa et al. 2019), and the clump epoch merger hypothesis, since the merger histories of galaxies are
lasts only a brief time, until the gas mass fraction declines very variable (e.g., Lacey & Cole 1993; Stewart et al. 2008;
(Cacciato et al. 2012). We showed in Paper I that the clumps that Boylan-Kolchin et al. 2010).
form in the clumpy simulation are comparable to the ones found Upcoming data from the James Webb Space Telescope will
in high-redshift galaxies and predicted that chemical bimodal- measure the chemistry of the Andromeda galaxy’s disk from
ities in disks should be common. Using MUSE spectroscopy, resolved stellar spectroscopy. Andromeda is known to have had
Scott et al. (2021) showed that the MW analog UGC 10738 has a much more active merger history than that of the MW (e.g.,
an α-rich geometric thick disk, from which they conclude that McConnachie et al. 2010; Weisz et al. 2014; D’Souza &
accretion events are unlikely sources of thick disks. More studies Bell 2018; Hammer et al. 2018; McConnachie et al. 2018). If
such as this can help establish whether geometric thick disks are clumps played an important role in its chemical evolution, we
α-enhanced. This will be particularly important for exploring the expect the chemistry of Andromeda’s old disks to be comprised
13
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
of, at least, a high-α and a low-α track somewhat similar to the by funnelling gas to the center, leading to further star formation
MW’s, with possibly additional merger-induced tracks. and compaction (Dekel & Burkert 2014).
However more detailed tests must necessarily come from the However other studies have questioned the importance of
MW since we can study it in much greater detail than any other clumps for the evolution of galaxies. Efficient coupling of
galaxy. Understanding to what extent the outcome of the GSE feedback energy to gas destroys clumps (Elmegreen et al. 2008;
merger is degenerate with the clump scenario is an important Hopkins et al. 2012) and many simulations that employ high
ingredient in unravelling the formation of the MW. A holistic feedback prescriptions have failed to find significant clumps or
approach, considering the properties of the bulge, the thin have found ones that do not contribute much to bulges (e.g.,
+thick disks, and the splash, is vital to this enterprise. However Tamburello et al. 2015). The short-lived clumps in the FIRE
efforts thus far have been hampered by the relatively small data simulations do not manage to migrate to the bulge (Oklopčić
sets comprising thousands of stars. Future space-based (e.g., et al. 2017), and may not even have been bound. Similarly, in
Gaia) and ground-based observatories’ (e.g., Vera Rubin the NIHAO simulation suite, Buck et al. (2017) found that
Telescope) surveys will permit proper-motion measurements clumps are only present in the light, not in the mass, and
of large samples of bulge stars to help unravel the formation of therefore have minimal contribution to bulge growth.
the bulge (Gough-Kelly et al. 2022). The detailed treatment of various forms of feedback (e.g.,
A possible test is the distribution of ages at the bulge’s low- Fensch & Bournaud 2021) therefore plays an important role in
[Fe/H] peak versus that of the trough between the two peaks. the ease with which clumps form in simulations. Moreover, in
Most stars in the MW’s bulge will now be old; measuring an simulations of single giant molecular clouds (GMCs), the
age difference of ∼2 Gyr in a present-day ∼10 Gyr old bulge energy imparted by photoionization, winds, and supernova
(Ortolani et al. 1995; Kuijken & Rich 2002; Ferreras et al. feedback can be channeled along preferred directions, thereby
2003; Zoccali et al. 2003; Sahu et al. 2006; Clarkson et al. preserving the GMC for a longer time than would otherwise be
2008, 2011; Brown et al. 2010; Valenti et al. 2013; Calamida expected (Rogers & Pittard 2013; Dale 2017; Howard et al.
et al. 2014; Renzini et al. 2018; Surot et al. 2019) is 2017). Thus tests of what role, if any, clumps have played in
challenging. Nonetheless, the chemical thin and thick disks the evolution of galaxies like the MW can inform improve-
do appear to overlap in age, as seen by the existence of ments in subgrid implementations of feedback on the smallest
RR Lyrae with small vertical excursions and low [α/Fe] scales, perhaps by accounting for feedback channeling. Further
(Prudil et al. 2020). An age overlap between the MW’s thin and study of the impact of clumps on galaxy formation therefore
thick disk, which was predicted in Paper I, has also been may have much broader impact on the study of galaxy
demonstrated by Beraldo e Silva et al. (2021) using the stellar formation (e.g., Dekel et al. 2022). Recently, Marasco et al.
ages of turnoff and giant stars from the Sanders & Das (2018) (2023) found that the observed outflows from a sample of
catalog. Likewise, Silva Aguirre et al. (2018) find an age starbursting dwarf galaxies are lower than predicted by
overlap between high-α and low-α disk stars from astroseismic cosmological simulations that employ high feedback. They
ages. Gent et al. (2022) reach a similar conclusion based on find mass loading factors of warm gas outflows more than 2
data from the Gaia-ESO survey together with Gaia EDR3 data. orders of magnitude lower than predicted, providing strong
Thus it may well be possible to measure the mean age support for the need for gentler feedback prescriptions.
difference between stars at the trough and those in the low-[Fe/
H] peak to test whether clumps have contributed to the bulge. 4.4. Caveats
The simulation presented in this paper is clearly idealized
4.3. Clumps as Probes of Feedback Implementations and lacks some of the ingredients that have been suggested to
have mattered in the MW’s chemical evolution. Of these the
Clumps were first proposed to play a role in the formation of most important is the merger of the GSE progenitor. The effect
bulges by Noguchi (1999). Following this suggestion, several of the GSE merger will be explored in future papers (e.g.,
works explored the role of clumps in bulge formation (Immeli Amarante et al. 2022).
et al. 2004; Bournaud et al. 2007; Elmegreen et al. 2008; Our simulations place the initial gas in a hot corona. This is
Aumer et al. 2010; Inoue & Saitoh 2012). When the appropriate for a galaxy of the MW’s mass since redshift z∼ 1
cosmological setting is also included, the possibility of (Birnboim & Dekel 2003; Kereš et al. 2005), but is less
“ex situ” clumps forming directly in the cold gas streaming appropriate before then. More realistically, the gas should flow
in before reaching the disk was also recognized (Dekel et al. in along filaments (cold-mode accretion). Unfortunately setting
2009; Ceverino et al. 2010). Clumps can even be excited by up such initial conditions for high-resolution simulations is
external perturbations (Inoue et al. 2016). The cosmological difficult. However cosmological simulations have found
simulations of Dubois et al. (2021) find that ∼10% of the stellar clumps in galaxies still in the filamentary cold accretion mode
mass of z= 4 galaxies may be in the form of clumps, while (Dekel et al. 2009; Ceverino et al. 2010); provided that the
those of Mandelker et al. (2014) resulted in 60% of galaxies inflow rate of the gas, and the resulting clump and star
hosting an in situ clump population. Meanwhile Mandelker formation rates are realistic, it does not matter how gas reaches
et al. (2017) showed that bulges can host their own clump, the disk. If gas stalls in the outer disk (for instance as some of
which is more robust to feedback; they further showed that the gas does in Agertz et al. 2021), then it may be that gas
radiation pressure increases the cold gas fraction (by delaying surface densities are never high enough for clumps to form. We
star formation), increasing the lifetime of low-mass clumps. speculate that if this inhibits the flow of gas to the bulge then
Inoue & Saitoh (2012) showed that the bulges formed from the bulge itself may never reach the same high-[α/Fe] state as
clump mergers are rapidly rotating, exponential, and comprised the thick disk. However in general filamentary, cold-mode
of old, metal-rich stars, similar to the bulge of the MW. In accretion need not alter the general picture much so long as
addition, clumps may further affect the formation of the bulge realistic clumps form.
14
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
One limitation of the clumpy model presented here is that the rich, α-poor stars that form in situ after the epoch of
clump population in this particular simulation may be too large. clump formation. Such twin peaks are not present when
This is suggested by high rotation velocity at the center (Clarke the feedback suppresses clump formation. The relative
et al. 2019), and the failure to form a bar (although bars often mass in the high- and low-[Fe/H] peaks constrains the
fail to form at this mass resolution). These effects may be epoch when star formation in the bulge is quenched (see
improved in models with higher feedback that still permit long- Sections 3.1, 3.2, and 3.3).
lived, but lower-mass clumps, which are still able to produce a 2. Star formation within the bulge occurs in the high-ΣSFR
high-α population (e.g., Garver et al. 2023). clump mode. This ensures that a separate low-[α/Fe]
Despite the absence of a bar, we may anticipate what the track never forms (see Section 3.4).
influence of a bar might be. Debattista et al. (2017) showed that 3. The metal-rich bulge population, while on average
many of the trends with metallicity observed in the MW’s younger than the metal-poor population, overlaps with
bulge can be explained by the secular evolution of the bar, via a it in age because the latter population is partly built from
mechanism they termed kinematic fractionation. In this stars that came in with clumps after the chemical
mechanism, different populations are separated by the bar evolution of in situ star formation in the bulge had
formation on the basis of their radial random motion. moved to higher metallicities (see Section 3.5).
Populations with large radial random motions are lifted by 4. By the end of the clump epoch, the bulge is already
the bar to large heights ending as a spheroidal population and rapidly rotating. The high-[α/Fe], low-[Fe/H] bulge
forming a weaker bar, whereas populations that start with lower population is kinematically hotter than the low-[α/Fe],
radial random motion do not rise to as large heights but end high-[Fe/H] one (see Section 3.6).
with a more strongly peanut-shaped distribution and a stronger 5. The population at the metal-rich peak is prominent at low
bar. As a result, in general the X-shape is better traced by the latitudes but declines with distance from the mid-plane,
metal-rich stars, which are younger and start out cooler, while as observed in the MW (see Section 3.7).
metal-poor stars, which are older and thus kinematically hotter, 6. A test of the role of clumps on the MW’s bulge comes
trace a more boxy structure (Athanassoula et al. 2017; from a comparison of the age distributions of the low-
Debattista et al. 2017; Buck et al. 2018; Debattista et al. [Fe/H] peak and the trough between the two peaks. In the
2019; Fragkoudi et al. 2020). Subsequently, Debattista et al. presence of clumps, the age distributions overlap
(2020) demonstrated that the vertical thickening of stellar significantly, with the mean age higher in the trough
populations increases monotonically with the radial action of than at the low-[Fe/H] peak, contrary to the usual
stars from before the bar formed. Since stars with larger radial expectation of increasing metallicity with age (see
random motion are typically older, and usually more metal- Section 3.5).
poor, kinematic fractionation results in a vertical metallicity
gradient. In addition the X-shape ends up better traced by This paper, together with Paper I, presents idealized
metal-rich stars, as is observed in the MW (Ness et al. 2012; simulations that demonstrate that clump formation provides a
Uttenthaler et al. 2012; Rojas-Arriagada et al. 2014). We have very direct and natural way of producing chemical trends
shown here that the radial action, J , decreases along the observed not only in the MW’s thin+thick disk but also in the
R
bulge’s chemical track with increasing metallicity. Thus bulge. The simulations are by no means wholly realistic, but
kinematic fractionation would raise stars at the metal-poor the ease with which they produce the trends observed in the
peak to larger heights than those of the metal-rich peak. This MW encourages us to explore further the role of clumps in the
would further enhance the trends of Section 3.7, which already early history of galaxies. In contrast, satisfying both constraints
match those observed in the MW. in other scenarios of thick disk formation may require a more
Recently Queiroz et al. (2021) derived distances of a large specific set of circumstances, which would mean the MW is
sample of bulge stars with STARHORSE using APOGEE DR16 unusual. The clump model makes some important predictions
and Gaia EDR3 parallaxes. They argued that the chemistry of that can be verified with future facilities, including that
the bulge is comprised of not one but two tracks, contrary to chemical thick disks should be common in MW-mass galaxies
earlier studies. The two tracks do not overlap in metallicity and that a population of chemical thin-disk stars of comparable
(unlike the thin+thick disks), but are separated by a gap and age to the thick disk should exist in the MW. Studying the
have different slopes. After accounting for the stellar popula- consequences of the clumps in such simulations may also
tion-dependent selection function of APOGEE, Eilers et al. provide a useful probe of feedback implementations.
(2022) also found tracks with different slope, although the
tracks still do not overlap in metallicity. If these trends are V.P.D., L.B.S., and T.K. were supported by STFC
confirmed by imminent large surveys with instruments such as Consolidated grant ST/R000786/1. D.J.L. was supported for
MOONS (Cirasuolo et al. 2014) and 4MOST (de Jong et al. part of this project by a UCLan UURIP internship. L.B.S
2014), this may suggest that the clump scenario needs acknowledges the support of NASA-ATP award
alteration or is perhaps wrong. 80NSSC20K0509 and National Science Foundation Astron-
omy and Astrophysics Research Grants (AAG) grant AST-
2009122. J.A.S.A. acknowledges funding from the European
4.5. Summary Research Council (ERC) under the European Unionʼs Horizon
The principal results of this paper are as follows: 2020 research and innovation program (grant agreement No.
852839). M.Z. acknowledges support from the ANID BASAL
1. A single track with two peaks in the bulge’s [Fe/H]–[α/ Center for Astrophysics and Associated Technologies (CATA)
Fe] space results when clump formation can occur in the through grants AFB170002, ACE210002, and FB210003, the
early evolution. Clumps sink to the center, contributing to ANID Millennium Institute of Astrophysics (MAS)
the bulge. The bulge is later populated by more metal- ICN12_009, and ANID Fondecyt Regular grant 1191505. E.
15
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
V. acknowledges the Excellence Cluster ORIGINS funded by Beraldo e Silva, L., Debattista, V. P., Khachaturyants, T., & Nidever, D. 2020,
the Deutsche Forschungsgemeinschaft (DFG; German MNRAS, 492, 4716
Research Foundation) under Germanyʼs Excellence Strategy Beraldo e Silva, L., Debattista, V. P., Nidever, D., & Amarante, J. A. S. 2021,
MNRAS, 502, 260
EXC-2094-390783311. S.A.K. would like to acknowledge Binney, J. 2012, MNRAS, 426, 1324
support from NASAʼs Astrophysics Data Analysis Program Bird, J. C., Kazantzidis, S., Weinberg, D. H., et al. 2013, ApJ, 773, 43
(ADAP) grant number 80NSSC20K0760. We thank the Bird, J. C., Loebman, S. R., Weinberg, D. H., et al. 2021, MNRAS, 503, 1815
anonymous referee for comments that helped improve this Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349
paper. An important part of the methodology for the stellar Bland-Hawthorn, J., & Gerhard, O. 2016, ARA&A, 54, 529
Bournaud, F., Elmegreen, B. G., & Elmegreen, D. M. 2007, ApJ, 670, 237
population modeling used in this paper was worked out in 2018 Bournaud, F., Elmegreen, B. G., & Martig, M. 2009, ApJL, 707, L1
at the Aspen Center for Physics, which is supported by Bournaud, F., Perret, V., Renaud, F., et al. 2014, ApJ, 780, 57
National Science Foundation grant PHY-1607611. The visit of Bovy, J., Rix, H.-W., Liu, C., et al. 2012, ApJ, 753, 148
V.P.D. was partially supported by a grant from the Simons Bovy, J., Rix, H.-W., Schlafly, E. F., et al. 2016, ApJ, 823, 30
Boylan-Kolchin, M., Springel, V., White, S. D. M., & Jenkins, A. 2010,
Foundation. The simulations in this paper were run at the MNRAS, 406, 896
DiRAC Shared Memory Processing system at the University of Brook, C. B., Gibson, B. K., Martel, H., & Kawata, D. 2005, ApJ, 630, 298
Cambridge, operated by the COSMOS Project at the Depart- Brown, T. M., Sahu, K., Anderson, J., et al. 2010, ApJL, 725, L19
ment of Applied Mathematics and Theoretical Physics on Buck, T. 2020, MNRAS, 491, 5435
behalf of the STFC DiRAC High Performance Computing Buck, T., Maccio, A. V., Obreja, A., et al. 2017, MNRAS, 468, 3628
Buck, T., Ness, M. K., Maccio, A. V., Obreja, A., & Dutton, A. A. 2018, ApJ,
(HPC) Facility: www.dirac.ac.uk. This equipment was funded 861, 88
by theDepartment for Business, Innovation and Skills (BIS) Bullock, J. S., Dekel, A., Kolatt, T. S., et al. 2001, ApJ, 555, 240
National E-infrastructure capital grant ST/J005673/1, STFC Cacciato, M., Dekel, A., & Genel, S. 2012, MNRAS, 421, 818
capital grant ST/H008586/1, and STFC DiRAC Operations Calamida, A., Sahu, K. C., Anderson, J., et al. 2014, ApJ, 790, 164
Cava, A., Schaerer, D., Richard, J., et al. 2018, NatAs, 2, 76
grant ST/K00333X/1. DiRAC is part of the National Ceverino, D., Dekel, A., & Bournaud, F. 2010, MNRAS, 404, 2151
E-Infrastructure. This paper is dedicated to the memory of Chiappini, C. 2009, in IAU Symp. 254, The Galaxy Disk in Cosmological
George Lake, for whose support and inspiration V.P.D. is Context, ed. J. Andersen, B. m. Nordströara, & J. Bland-Hawthorn
deeply indebted. (Cambridge: Cambridge Univ. Press), 191
Chiappini, C., Matteucci, F., & Gratton, R. 1997, ApJ, 477, 765
Cirasuolo, M., Afonso, J., Carollo, M., et al. 2014, Proc. SPIE, 9147, 91470N
ORCID iDs Clarke, A. J., Debattista, V. P., Nidever, D. L., et al. 2019, MNRAS, 484, 3476
Clarkson, W., Sahu, K., Anderson, J., et al. 2008, ApJ, 684, 1110
Victor P. Debattista https://orcid.org/0000-0001-7902-0116 Clarkson, W. I., Sahu, K. C., Anderson, J., et al. 2011, ApJ, 735, 37
Oscar A. Gonzalez https://orcid.org/0000-0003-2478-6020 Conroy, C., Weinberg, D. H., Naidu, R. P., et al. 2022, arXiv:2204.02989
Leandro Beraldo e Silva https://orcid.org/0000-0002- Dale, J. E. 2017, MNRAS, 467, 1067
0740-1507 de Jong, R. S., Barden, S., Bellido-Tirado, O., et al. 2014, Proc. SPIE, 9147,
91470M
João A. S. Amarante https://orcid.org/0000-0002- Debattista, V. P., Gonzalez, O. A., Sanderson, R. E., et al. 2019, MNRAS,
7662-5475 485, 5073
Manuela Zoccali https://orcid.org/0000-0002-5829-2267 Debattista, V. P., Liddicott, D. J., Khachaturyants, T., & Beraldo e Silva, L.
Elena Valenti https://orcid.org/0000-0002-6092-7145 2020, MNRAS, 498, 3334
Debattista, V. P., Ness, M., Gonzalez, O. A., et al. 2017, MNRAS, 469, 1587
Deanne B. Fisher https://orcid.org/0000-0003-0645-5260 Dekel, A., Birnboim, Y., Engel, G., et al. 2009, Natur, 457, 451
Tigran Khachaturyants https://orcid.org/0000-0002- Dekel, A., & Burkert, A. 2014, MNRAS, 438, 1870
3343-6615 Dekel, A., Mandelker, N., Bournaud, F., et al. 2022, MNRAS, 511, 316
David L. Nidever https://orcid.org/0000-0002-1793-3689 Dessauges-Zavadsky, M., Schaerer, D., Cava, A., Mayer, L., & Tamburello, V.
Thomas R. Quinn https://orcid.org/0000-0001-5510-2803 2017, ApJL, 836, L22
Di Matteo, P. 2016, PASA, 33, 027
Min Du https://orcid.org/0000-0001-9953-0359 Di Matteo, P., Fragkoudi, F., Khoperskov, S., et al. 2019a, A&A, 628, A11
Susan Kassin https://orcid.org/0000-0002-3838-8093 Di Matteo, P., Haywood, M., Lehnert, M. D., et al. 2019b, A&A, 632, A4
D’Souza, R., & Bell, E. F. 2018, NatAs, 2, 737
References Dubois, Y., Beckmann, R., Bournaud, F., et al. 2021, A&A, 651, A109
Eilers, A.-C., Hogg, D. W., Rix, H.-W., et al. 2022, ApJ, 928, 23
Elmegreen, B. G., Bournaud, F., & Elmegreen, D. M. 2008, ApJ, 688, 67
Agertz, O., Renaud, F., Feltzing, S., et al. 2021, MNRAS, 503, 5826 Elmegreen, B. G., & Elmegreen, D. M. 2005, ApJ, 627, 632
Agertz, O., Teyssier, R., & Moore, B. 2009, MNRAS, 397, L64 Elmegreen, D. M., Elmegreen, B. G., Ravindranath, S., & Coe, D. A. 2007,
Alves-Brito, A., Melendez, J., Asplund, M., Ramirez, I., & Yong, D. 2010, ApJ, 658, 763
A&A, 513, A35 Fensch, J., & Bournaud, F. 2021, MNRAS, 505, 3579
Amarante, J. A. S., Beraldo e Silva, L., Debattista, V. P., & Smith, M. C. 2020, Ferreras, I., Wyse, R. F. G., & Silk, J. 2003, MNRAS, 345, 1381
ApJL, 891, L30 Fisher, D. B., Glazebrook, K., Abraham, R. G., et al. 2017, ApJL, 839, L5
Amarante, J. A. S., Debattista, V. P., Beraldo E Silva, L., Laporte, C. F. P., & Forbes, D. A. 2020, MNRAS, 493, 847
Deg, N. 2022, ApJ, 937, 12 Forster Schreiber, N. M., Shapley, A. E., Genzel, R., et al. 2011, ApJ, 739, 45
Ambachew, L., Fisher, D. B., Glazebrook, K., et al. 2022, MNRAS, 512, 3079
Fragkoudi, F., Grand, R. J. J., Pakmor, R., et al. 2020, MNRAS, 494, 5936
Athanassoula, E., Rodionov, S. A., & Prantzos, N. 2017, MNRAS, 467, L46
Aumer, M., Burkert, A., Johansson, P. H., & Genzel, R. 2010, ApJ, 719, 1230 Freeman, K., Ness, M., Wylie-de-Boer, E., et al. 2013, MNRAS, 428, 3660
Bekki, K., & Tsujimoto, T. 2011, MNRAS, 416, L60 Gallart, C., Bernard, E. J., Brook, C. B., et al. 2019, NatAs, 3, 932
Belokurov, V., Erkal, D., Evans, N. W., Koposov, S. E., & Deason, A. J. 2018, Garver, B. R., Nidever, D. L., Debattista, V. P., Beraldo e Silva, L., &
MNRAS, 478, 611 Khachaturyants, T. 2023, ApJ, submitted
Belokurov, V., & Kravtsov, A. 2022, MNRAS, 514, 689 Genel, S., Naab, T., Genzel, R., et al. 2012, ApJ, 745, 11
Belokurov, V., Sanders, J. L., Fattahi, A., et al. 2020, MNRAS, 494, 3880 Gent, M. R., Eitner, P., Laporte, C. F. P., et al. 2022, arXiv:2206.10949
Benincasa, S. M., Wadsley, J. W., Couchman, H. M. P., Pettitt, A. R., & Genzel, R., Newman, S., Jones, T., et al. 2011, ApJ, 733, 101
Tasker, E. J. 2019, MNRAS, 486, 5022 Gonzalez, O. A., Zoccali, M., Vasquez, S., et al. 2015, A&A, 584, A46
Bensby, T., Feltzing, S., Gould, A., et al. 2017, A&A, 605, A89 Gough-Kelly, S., Debattista, V. P., Clarkson, W. I., et al. 2022, MNRAS,
Bensby, T., Feltzing, S., Johnson, J. A., et al. 2010, A&A, 512, A41 509, 4829
Bensby, T., Yee, J. C., Feltzing, S., et al. 2013, A&A, 549, A147 Grand, R. J. J., Bustamante, S., Gomez, F. A., et al. 2018, MNRAS, 474, 3629
16
The Astrophysical Journal, 946:118 (17pp), 2023 April 1 Debattista et al.
Grand, R. J. J., Gómez, F. A., Marinacci, F., et al. 2017, MNRAS, 467, 179 Ness, M., Freeman, K., Athanassoula, E., et al. 2012, ApJ, 756, 22
Grisoni, V., Spitoni, E., Matteucci, F., et al. 2017, MNRAS, 472, 3637 Ness, M., Freeman, K., Athanassoula, E., et al. 2013, MNRAS, 430, 836
Guo, Y., Ferguson, H. C., Bell, E. F., et al. 2015, ApJ, 800, 39 Nidever, D. L., Bovy, J., Bird, J. C., et al. 2014, ApJ, 796, 38
Guo, Y., Giavalisco, M., Ferguson, H. C., Cassata, P., & Koekemoer, A. M. Noguchi, M. 1999, ApJ, 514, 77
2012, ApJ, 757, 120 Oklopčić, A., Hopkins, P. F., Feldmann, R., et al. 2017, MNRAS, 465, 952
Guo, Y., Rafelski, M., Bell, E. F., et al. 2018, ApJ, 853, 108 Orkney, M. D. A., Laporte, C. F. P., Grand, R. J. J., et al. 2022, MNRAS,
Hammer, F., Yang, Y. B., Wang, J. L., et al. 2018, MNRAS, 475, 2754 517, L138
Hartmann, M., Debattista, V. P., Seth, A., Cappellari, M., & Quinn, T. R. 2011, Ortolani, S., Renzini, A., Gilmozzi, R., et al. 1995, Natur, 377, 701
MNRAS, 418, 2697 Prudil, Z., Dekany, I., Grebel, E. K., & Kunder, A. 2020, MNRAS, 492, 3408
Hayden, M. R., Bovy, J., Holtzman, J. A., et al. 2015, ApJ, 808, 132 Queiroz, A. B. A., Chiappini, C., Perez-Villegas, A., et al. 2021, A&A,
Helmi, A., Babusiaux, C., Koppelman, H. H., et al. 2018, Natur, 563, 85 656, A156
Hill, V., Lecureur, A., Gomez, A., et al. 2011, A&A, 534, A80 Ravindranath, S., Giavalisco, M., Ferguson, H. C., et al. 2006, ApJ, 652, 963
Hopkins, P. F., Keres, D., Murray, N., Quataert, E., & Hernquist, L. 2012, Read, J. I., Lake, G., Agertz, O., & Debattista, V. P. 2008, MNRAS, 389, 1041
MNRAS, 427, 968 Renaud, F., Agertz, O., Andersson, E. P., et al. 2021a, MNRAS, 503, 5868
Horta, D., Schiavon, R. P., Mackereth, J. T., et al. 2021, MNRAS, 500, 1385 Renaud, F., Agertz, O., Read, J. I., et al. 2021b, MNRAS, 503, 5846
Howard, C. S., Pudritz, R. E., & Harris, W. E. 2017, MNRAS, 470, 3346 Renzini, A., Gennaro, M., Zoccali, M., et al. 2018, ApJ, 863, 16
Huertas-Company, M., Guo, Y., Ginzburg, O., et al. 2020, MNRAS, 499, 814 Rich, R. M. 1988, AJ, 95, 828
Immeli, A., Samland, M., Gerhard, O., & Westera, P. 2004, A&A, 413, 547 Rix, H.-W., Chandra, V., Andrae, R., et al. 2022, ApJ, 941, 45
Inoue, S., Dekel, A., Mandelker, N., et al. 2016, MNRAS, 456, 2052 Rogers, H., & Pittard, J. M. 2013, MNRAS, 431, 1337
Inoue, S., & Saitoh, T. R. 2012, MNRAS, 422, 1902 Rojas-Arriagada, A., Recio-Blanco, A., Hill, V., et al. 2014, A&A, 569, A103
Inoue, S., & Yoshida, N. 2019, MNRAS, 488, 4400 Rojas-Arriagada, A., Zasowski, G., Schultheis, M., et al. 2020, MNRAS,
Johnson, C. I., Rich, R. M., Kobayashi, C., Kunder, A., & Koch, A. 2014, AJ, 499, 1037
148, 67 Rojas-Arriagada, A., Zoccali, M., Schultheis, M., et al. 2019, A&A, 626, A16
Johnson, C. I., Rich, R. M., Simion, I. T., et al. 2022, MNRAS, 515, 1469 Sahu, K. C., Casertano, S., Bond, H. E., et al. 2006, Natur, 443, 534
Jonsson, H., Ryde, N., Schultheis, M., & Zoccali, M. 2017, A&A, 598, A101 Sanders, J. L., & Das, P. 2018, MNRAS, 481, 4093
Kassin, S. A., Weiner, B. J., Faber, S. M., et al. 2012, ApJ, 758, 106 Schultheis, M., Rojas-Arriagada, A., Garcia Perez, A. E., et al. 2017, A&A,
Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2 600, A14
Khachaturyants, T., Beraldo e Silva, L., & Debattista, V. P. 2021, MNRAS, Scott, N., van de Sande, J., Sharma, S., et al. 2021, ApJL, 913, L11
508, 2350 Shen, J., Rich, R. M., Kormendy, J., et al. 2010a, ApJL, 720, L72
Khoperskov, S., Haywood, M., Di Matteo, P., Lehnert, M. D., & Combes, F. Shen, S., Wadsley, J., & Stinson, G. 2010b, MNRAS, 407, 1581
2018, A&A, 609, A60 Silva Aguirre, V., Bojsen-Hansen, M., Slumstrup, D., et al. 2018, MNRAS,
Khoperskov, S., Haywood, M., Snaith, O., et al. 2021, MNRAS, 501, 5176 475, 5487
Kruijssen, J. M. D., Pfeffer, J. L., Chevance, M., et al. 2020, MNRAS, Silverman, B. W. 1986, Density Estimation for Statistics and Data Analysis
498, 2472 (London: Chapman and Hall)
Kuijken, K., & Rich, R. M. 2002, AJ, 124, 2054 Snaith, O. N., Bailin, J., Gibson, B. K., et al. 2016, MNRAS, 456, 3119
Lacey, C., & Cole, S. 1993, MNRAS, 262, 627 Spitoni, E., Verma, K., Silva Aguirre, V., et al. 2021, A&A, 647, A73
Lane, J. M. M., Bovy, J., & Mackereth, J. T. 2022, MNRAS, 510, 5119 Stewart, K. R., Bullock, J. S., Wechsler, R. H., Maller, A. H., & Zentner, A. R.
Lenkic, L., Bolatto, A. D., Fisher, D. B., et al. 2021, MNRAS, 506, 3916 2008, ApJ, 683, 597
Lian, J., Zasowski, G., Hasselquist, S., et al. 2020, MNRAS, 497, 3557 Stinson, G., Seth, A., Katz, N., et al. 2006, MNRAS, 373, 1074
Lian, J., Zasowski, G., Mackereth, T., et al. 2022, MNRAS, 513, 4130 Surot, F., Valenti, E., Hidalgo, S. L., et al. 2019, A&A, 623, A168
Livermore, R. C., Jones, T., Richard, J., et al. 2012, MNRAS, 427, 688 Tamburello, V., Mayer, L., Shen, S., & Wadsley, J. 2015, MNRAS, 453, 2490
Livermore, R. C., Jones, T. A., Richard, J., et al. 2015, MNRAS, 450, 1812 Tamfal, T., Mayer, L., Quinn, T. R., et al. 2022, ApJ, 928, 106
Mackereth, J. T., Crain, R. A., Schiavon, R. P., et al. 2018, MNRAS, 477, 5072 Ted Mackereth, J., Bovy, J., Schiavon, R. P., et al. 2017, MNRAS, 471, 3057
Mackereth, J. T., Schiavon, R. P., Pfeffer, J., et al. 2019, MNRAS, 482, 3426 Tsujimoto, T., & Bekki, K. 2012, ApJ, 747, 125
Majewski, S. R., Schiavon, R. P., Frinchaboy, P. M., et al. 2016, AN, 337, 863 Uttenthaler, S., Schultheis, M., Nataf, D. M., et al. 2012, A&A, 546, A57
Mandelker, N., Dekel, A., Ceverino, D., et al. 2014, MNRAS, 443, 3675 Valenti, E., Zoccali, M., Renzini, A., et al. 2013, A&A, 559, A98
Mandelker, N., Dekel, A., Ceverino, D., et al. 2017, MNRAS, 464, 635 Vasiliev, E. 2019, MNRAS, 482, 1525
Marasco, A., Belfiore, F., Cresci, G., et al. 2023, A&A, 670, A92 Vincenzo, F., & Kobayashi, C. 2020, MNRAS, 496, 80
Massari, D., Koppelman, H. H., & Helmi, A. 2019, A&A, 630, L4
Wadsley, J. W., Stadel, J., & Quinn, T. 2004, NewA, 9, 137
McConnachie, A. W., Ferguson, A. M. N., Irwin, M. J., et al. 2010, ApJ,
723, 1038 Weisz, D. R., Skillman, E. D., Hidalgo, S. L., et al. 2014, ApJ, 789, 24
McConnachie, A. W., Ibata, R., Martin, N., et al. 2018, ApJ, 868, 55 Williams, A. A., Evans, N. W., Molloy, M., et al. 2016, ApJL, 824, L29
McWilliam, A., & Rich, R. M. 1994, ApJS, 91, 749 Wisnioski, E., Forster Schreiber, N. M., Wuyts, S., et al. 2015, ApJ, 799, 209
Melendez, J., Asplund, M., Alves-Brito, A., et al. 2008, A&A, 484, L21 Wuyts, S., Forster Schreiber, N. M., Genzel, R., et al. 2012, ApJ, 753, 114
Meng, X., & Gnedin, O. Y. 2021, MNRAS, 502, 1433 Wylie, S. M., Gerhard, O. E., Ness, M. K., et al. 2021, A&A, 653, A143
Myeong, G. C., Belokurov, V., Aguado, D. S., et al. 2022, ApJ, 938, 21 Zasowski, G., Schultheis, M., Hasselquist, S., et al. 2019, ApJ, 870, 138
Naidu, R. P., Ji, A. P., Conroy, C., et al. 2022, ApJL, 926, L36 Zoccali, M., Gonzalez, O. A., Vasquez, S., et al. 2014, A&A, 562, A66
Nataf, D. M. 2017, PASA, 34, e041 Zoccali, M., Renzini, A., Ortolani, S., et al. 2003, A&A, 399, 931
Navarro, J. F., Frenk, C. S., & White, S. D. M. 1997, ApJ, 490, 493 Zoccali, M., Vasquez, S., Gonzalez, O. A., et al. 2017, A&A, 599, A12
17