# The proto-galaxy of Milky Way-mass haloes in the FIRE simulations

Danny Horta1, Emily C. Cunningham2,1†, Robyn Sanderson1,3, Kathryn V.
Johnston2, Alis Deason4,5, Andrew Wetzel6, Fiona McCluskey6, Nicolás Garavito-
Camargo1, Lina Necib7,8, Claude-André Faucher-Giguère9, Arpit Arora3, Pratik
J. Gandhi6
1Center for Computational Astrophysics, Flatiron Institute, 162 5th Ave., New
York, NY 10010, USA
2Department of Astronomy, Columbia University, 550 West 120th Street, New
York, NY 10027, USA
3Department of Physics $\&$ Astronomy, University of Pennsylvania, 209 S 33rd
St., Philadelphia, PA 19104, USA
4Institute for Computational Cosmology, Department of Physics, Durham
University, South Road, Durham DH1 3LE, UK
5Centre for Extragalactic Astronomy, Department of Physics, Durham University,
South Road, Durham DH1 3LE, UK
6Department of Physics and Astronomy, University of California, Davis, CA
95616, USA
7The NSF AI Institute for Artificial Intelligence and Fundamental Interactions
8Department of Physics and Kavli Institute for Astrophysics and Space
Research, Massachusetts Institute of Technology,
77 Massachusetts Ave, Cambridge MA 02139, USA
9Department of Physics $\&$ Astronomy and CIERA, Northwestern University, 1800
Sherman Ave, Evanston, IL 60201, USA
† Hubble Fellow E-mail<EMAIL_ADDRESS>
(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

Observational studies are finding stars believed to be relics of the earliest
stages of hierarchical mass assembly of the Galaxy. In this work, we
contextualize these findings by studying the masses, ages, spatial
distributions, morphology, kinematics, and chemical compositions of proto-
galaxy populations from the 13 Milky Way (MW)-mass galaxies from the FIRE-2
cosmological zoom-in simulations. Our findings indicate that proto-Milky Way
populations: i) are predominantly centrally concentrated, with $\sim 50\%$ of
the stars contained within 5–10 kpc; ii) on average show weak but systematic
net rotation in the plane of the host’s disc at $z=0$ (i.e.,
0.25$\lesssim$$\langle$$\kappa$/$\kappa_{\mathrm{disc}}\rangle\lesssim 0.8$);
iii) present [$\alpha$/Fe]-[Fe/H] compositions that overlap with the metal-
poor tail of the host’s old disc; iv) tend to assemble slightly earlier in
Local Group-like environments than in systems in isolation. Interestingly, we
find that $\sim$60$\%$ of the proto-Milky Way galaxies are comprised by 1
dominant system (1/5$\lesssim$M⋆/M⋆,proto-MilkyWay$\lesssim$4/5) and 4–5 lower
mass systems (M⋆/M⋆,proto-MilkyWay$\lesssim$1/10); the other $\sim$40$\%$ are
comprised by 2 dominant systems and 3–4 lower mass systems. These
massive/dominant proto-Milky Way fragments can be distinguished from the lower
mass ones in chemical-kinematic samples, but appear (qualitatively)
indistinguishable from one another. Our results suggest that large/rich
chemical-kinematic-age samples of metal-poor stars in the inner Galaxy should
help characterise the different mass fragments of the proto-Milky Way. These
data may also help reveal if the Milky Way formed from one or two dominant
systems.

###### keywords:

Galaxy: general; Galaxy: formation; Galaxy: evolution; Galaxy: halo; Galaxy:
abundances; Galaxy: kinematics and dynamics

††pubyear: 2015††pagerange: The proto-galaxy of Milky Way-mass haloes in the
FIRE simulations–The proto-galaxy of Milky Way-mass haloes in the FIRE
simulations

##  1 Introduction

In the current accepted model for cosmology ($\Lambda$CDM), haloes grow by
accumulating lower mass building blocks through a process commonly referred to
as hierarchical mass assembly (White & Rees, 1978; White & Frenk, 1991,
e.g.,). During this process, the baryonic components (i.e., gas and stars that
constitute observable galaxies) of haloes also undergo this mechanism of mass
accumulation. This procedure is ubiquitous across the Universe and affects all
galaxies. Thus, the accretion history of a galaxy is a pivotal dictating
factor for its evolution and assembly of mass over time, and the clues for
disentangling such intricate process are all contained in the stellar halo.

Our current picture for the formation of stellar haloes suggests that they
form via a dual process. On the one hand, gas accretion from cosmic filaments
drives secular evolution and in situ star formation. On the other, the
accretion of different mass building blocks111A definition of building block
and main branch systems is provided in Table 1., each of which donate their
gas and stars to the resulting larger mass host, contribute in mass to a given
galaxy after becoming consumed. For the case of the Milky Way, early
observational evidence has suggested that the former process dominates within
the inner regions (i.e., $r$$\lesssim 20-30$ kpc), whereas the latter
dominates the outer regime (Chiba & Beers, 2000; Carollo et al., 2007; Deason
et al., 2011, e.g.,). However, more recent observational results suggest that
the in situ component may be a heated primordial disc (Bonaca et al., 2017; Di
Matteo et al., 2019; Belokurov et al., 2020, e.g.,). This dual formation
channel has also been somewhat shown to be the case in Andromeda (M31;
Ferguson et al., 2002; Brown et al., 2006; Escala et al., 2019, 2020, e.g.,).
The advent of detailed semi-analytic cosmological models and more recent
detailed cosmological simulations have supported this hypothesis on a
theoretical basis (Bullock & Johnston, 2005; Abadi et al., 2006; Bell et al.,
2008; Font et al., 2008; Johnston et al., 2008; Cooper et al., 2010; Font et
al., 2011; McCarthy et al., 2012; Amorisco, 2017; Khoperskov et al., 2022a,
b). However, whilst this picture may appear clear at current time, detailed
spectroscopic, photometric, and astrometric observations of halo stars with
Gaia (Gaia Collaboration et al., 2022) and large stellar surveys are starting
to poke holes in this framework.

Recent observational results have shown that the inner $\sim 30$ kpc of the
Milky Way’s stellar halo is awash with debris from engulfed satellite systems.
More specifically, it has been shown that in the Milky Way: $i$) the local
stellar halo is dominated by the debris of a massive and (likely) ancient
merged galaxy (Gaia-Sausage/Enceladus; Belokurov et al., 2018; Helmi et al.,
2018; Haywood et al., 2018; Mackereth et al., 2019), although see a counter-
argument by Donlon et al. (2022); $ii$) the innermost regions of the Galaxy
($r<5$ kpc) likely host the debris from an ancient and massive building block
(Heracles; Horta et al., 2021a)222It has been postulated that globular
clusters in the inner Galaxy could be linked to this system (Kruijssen et al.,
2020; Forbes, 2020).; this system could be related to another population
recently identified, referred to as the "proto-Milky Way" or Aurora (Belokurov
& Kravtsov, 2022; Conroy et al., 2021; Rix et al., 2022); $iii$) there are
numerous smaller-mass halo substructures/streams postulated to be the debris
from smaller mass accreted systems (see Helmi, 2020 for a review). Along those
lines, further out in the stellar halo, there are also streams from
currently/recently disrupted satellites (e.g., Sagittarius dSph, Ibata et al.,
1994; Cetus, Newberg et al., 2009; Orphan-Chenab, Grillmair & Dionatos, 2006;
Belokurov et al., 2007) 333And of course the Magellanic Clouds.. All these
different systems suggest that many components contribute to the formation of
the Milky Way’s stellar halo at all radii, in agreement with recent findings
using cosmological simulations (Khoperskov et al., 2022b; Horta et al., 2023b;
Orkney et al., 2023, e.g.,).

The recent plethora of observational findings hint towards the epoch of
formation for the Milky Way being early ($\gtrsim 8$ Gyr). This implies that
the major building blocks of the Galaxy are many gigayears old. While lower-
mass/unmixed accretions are interesting for unveiling the recent mass assembly
history of the Galaxy, or are useful in terms of potential measuring and
dynamical studies, it is the early building blocks and accretion events that:
$i$) constitute the bulk of the inner regions of the present day stellar halo
(i.e., $r\lesssim 30$ kpc); $ii$) supply the gas and baryonic material to fuel
the formation of the disc; $iii$) dictate key episodes in the formation of the
Milky Way. For this reason, observational studies have set out to investigate
the properties of the oldest stellar populations in the Milky Way, and have
reported chemical-kinematic evidence for the presence of stellar populations
in the innermost regions of the Galaxy –where one would expect the oldest
stars formed in situ to inhabit (El-Badry et al., 2018; Fragkoudi et al.,
2020)– that are distinguishable from the dominant bar/disc, and are likely to
constitute the entirety, or part of, the “proto-Milky Way”. This stellar
population has been postulated to arise from the main progenitor system of the
Milky Way, a major building block, or both simultaneously. Whilst the nature
of this population is yet to be fully established, these findings have
instigated the search for answers to additional pivotal questions when it
comes to understanding the formation of the Galaxy:

  * •

When did the proto-Milky Way form?

  * •

What constitutes the proto-Milky Way?

  * •

Is it possible to distinguish the different systems that formed the proto-
Milky Way from one another?

Santistevan et al. (2020) set out to investigate the formation times and
building blocks of Milky Way-mass galaxies in the FIRE-2 simulations, which is
directly related to the formation of proto-Milky Way populations. In their
work, the authors traced the origin of star particles formed at all times in
different regions of the Milky Way-mass galaxies, and found that: $i$) Milky
Way-mass galaxies typically assemble at $z\sim$3–4 (i.e., 11.6–12.2 Gyr ago);
$ii$) Milky Way-mass galaxies at $z=0$ are formed typically from $\sim$100
building blocks with M${}_{\star}\gtrsim 10^{5}$ M⊙; Milky Way-mass galaxies
in Local Group environments typically assemble earlier than those in
isolation, highlighting the importance of environment.

In this article, we set out to answer when and how proto-Milky Way systems may
have assembled, and what are the (present day) observable properties of
stellar populations that comprise a proto-Milky Way. To do so, we employ the
thirteen Milky Way-mass galaxies from the Latte and ELVIS suites (Wetzel et
al., 2016) of high-resolution cosmological simulations from the Feedback In
Realistic Environments (FIRE: Hopkins et al., 2018b) project (Section 2). We
take these simulations and track the star particles belonging to all the
resolvable luminous subhaloes, in order to identify the star particles
belonging to all the galactic systems that coalesce to form a proto-Milky Way
(Section 3). Upon selecting these populations, in Section 4 we go on to
examine the present day observable properties of their stellar debris (namely,
mass, age, spatial distribution, morphology, kinematics, and [$\alpha$/Fe] and
[Fe/H] chemical compositions), in order to test if different systems
constituting proto-Milky Way populations are distinguishable, and to provide
clues on how to identify proto-Milky Way populations. We then discuss our
results in the context of the current/future work in Section 5, provide the
limitations and improvements to our work in Section 6, and list our concluding
statements in Section 7.

##  2 Simulations

We make use of thirteen Milky Way-mass galaxies from the Latte (Wetzel et al.,
2016) and ELVIS (Garrison-Kimmel et al., 2019b) suites of FIRE-2444FIRE
project website: http://fire.northwestern.edu cosmological zoom-in
simulations. In detail, Latte is a suite of seven isolated Milky Way-mass
galaxies, and ELVIS is a suite of three Local Group-like pairs of Milky Way-
mass galaxies. Both these suites of simulations were run with the FIRE-2
physics model (Hopkins et al., 2018b), utilising the Lagrangian meshless
finite-mass $N$-body gravitational plus hydrodynamics code
GIZMO555http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html. (Hopkins,
2015). FIRE-2 simulations model many radiative cooling and heating processes
for gas, including free-free emission, photoionization/recombination, Compton
scattering, photoelectric, metal-line, molecular, fine-structure, dust-
collisional, and cosmic-ray heating across a temperature range of 10 - 1010K.
These simulations also include the spatially uniform, redshift-dependent,
cosmic UV background from Faucher-Giguère et al. (2009), for which HI
reionization occurs at $z_{\mathrm{reion}}$ $\sim$ 10\. Moreover, FIRE-2 self-
consistently generates and tracks 11 chemical abundance species (namely, H,
He, C, N, O, Ne, Mg, Si, S, Ca, and Fe), including sub-grid diffusion of these
abundances in gas via turbulence (Hopkins, 2016; Su et al., 2017; Escala et
al., 2018), as well as enrichment from core-collapse supernovae, supernovae
type Ia, and stellar winds.

The Latte suite has an initial baryonic mass resolution of 7,100 M⊙ for gas
particles, whereas the ELVIS suite has an initial baryonic mass resolution of
3,500-4000 M⊙. Latte uses a fixed resolution of 35,000 M⊙ for dark matter
particles, whereas ELVIS uses a fixed resolution of 19,000 M⊙. In both
Latte/ELVIS star formation occurs in gas that is self-gravitating, Jeans-
unstable, cold (T $<$ 104K), dense ($n$ $>$ 1,000 cm-3), and molecular
(following Krumholz & Gnedin 2011). Each star particle inherits the mass and
chemical abundance composition of its progenitor gas particle, and represents
a single stellar population with a Kroupa (2001) initial mass function. This
population evolves according to the STARBURST v7.0 models (Leitherer et al.,
1999), so the star particle decreases in mass with time as the most massive
stars die. By the present day most star particles have a mass of 4,000–5,000
M⊙. Stellar evolution gives rise to localized feedback at the location of each
star particle, including core-collapse and Ia supernovae, mass loss from
stellar winds, and radiation, including radiation pressure, photoionization,
and photo-electric heating. Implementation of these processes follows the
descriptions in Hopkins et al. (2018b) and Hopkins et al. (2018a).

The Latte/ELVIS suites were generated within periodic cosmological boxes of
lengths 70.4–172 Mpc using the code MUSIC (Hahn & Abel, 2011), and
constructing cosmological zoom-in initial conditions for each simulation at
$z$ $\simeq$ 99\. Each simulation has 600 snapshots saved down to $z$ = 0,
spaced every $\simeq$ 25 Myr. All simulations assume flat $\Lambda$-CDM
cosmology with parameters consistent with Planck Collaboration et al. (2020).
More specifically, the Latte suite (excluding m12r and m12w) used
$\Omega_{\mathrm{m}}$ = 0.272, $\Omega_{\mathrm{b}}$ = 0.0455,
$\sigma_{\mathrm{8}}$ = 0.807, $n_{\mathrm{s}}$ = 0.961, $h$ = 0.702. The m12r
and m12w halos were selected specifically because they host an LMC-mass
satellite galaxy, and they adopted more up-to-date initial conditions from
Planck Collaboration et al. (2020) compared to the rest of the Latte suite:
$h$ = 0.68, $\Omega_{\mathrm{\Lambda}}$ = 0.31, $\Omega_{\mathrm{m}}$ = 0.31,
$\Omega_{\mathrm{b}}$ = 0.048, $\sigma_{\mathrm{8}}$ = 0.82, $n_{\mathrm{s}}$
= 0.961. For the ELVIS suite, Thelma $\&$ Louise and Romulus $\&$ Remus both
used the same cosmology as in the original ELVIS dark matter only suite:
$\Omega_{\mathrm{m}}$ = 0.266, $\Omega_{\mathrm{b}}$ = 0.0449,
$\sigma_{\mathrm{8}}$ = 0.801, $n_{\mathrm{s}}$ = 0.963, $h$ = 0.71.
Conversely, Romeo $\&$ Juliet used the same cosmology as m12r/m12w. The post-
processing is done using gizmo analysis (Wetzel & Garrison-Kimmel, 2020b) and
halo analysis (Wetzel & Garrison-Kimmel, 2020a). Furthermore, dark matter
particles in each snapshot are processed with Rockstar (Behroozi et al., 2013)
to produce halo catalogs. For further details of how the Latte and ELVIS
suites were generated, we refer the reader to Wetzel et al. (2016) and
Garrison-Kimmel et al. (2019b) and references therein, respectively.

The resolution of this suite of simulations enables luminous subhalos to be
well resolved even near each Milky Way-like galaxy. It also resolves the
formation of tidal streams from satellite galaxies down to approximately 108M⊙
in total mass or 106M⊙ in stellar mass (at $z$=0), similar to that of the
Milky Way’s “classical" dwarf spheroidals (Panithanpaisal et al., 2021;
Cunningham et al., 2022; Shipp et al., 2022; Horta et al., 2023b, e.g.,).

The properties of the host galaxies in Latte show broad agreement with the
Milky Way, including the stellar-to-halo mass relation (Hopkins et al.,
2018b), stellar halos (Bonaca et al., 2017; Sanderson et al., 2018, 2020), and
the radial and vertical structure of their disks (Ma et al., 2017; Bellardini
et al., 2021). Moreover, the satellite populations of these simulation suites
have also been demonstrated to agree with several observed properties, such
as: the mass and velocity dispersions (Wetzel et al., 2016; Garrison-Kimmel et
al., 2019a); star formation histories (Garrison-Kimmel et al., 2019b); and
radial distributions (Samuel et al., 2020). Despite great similarities in the
properties of the satellite galaxies in Latte and that of the observed
satellites around the Milky Way, it has been shown that the former are
generally too metal-poor when compared to the latter (Escala et al., 2018;
Wheeler et al., 2019; Panithanpaisal et al., 2021). This could be in part
because of the assumed delay time distributions assumed for SNIa (Gandhi et
al., 2022), or because of the lack of modeling of Pop III stars. For this
work, we emphasize that we do not require quantitative agreement with
simulated and observed abundances. The relations found in this paper between
(luminous) haloes and their respective chemical abundances should be treated
qualitatively, and are intended for use within the simulations only. However,
although the _normalization_ of the various abundances is not always in good
agreement with observations, we expect the _trends_ we identify in this paper
to be robust.

##  3 Defining a proto-Milky Way

Table 1: List of definitions used in this article. See also Fig 1. Name | Definition
---|---
$t_{\mathrm{MR_{3:1}}}$ | Lookback time in the simulation we define the proto-Milky Way to emerge, defined as the time
| the main branch reaches a stellar mass ratio of 3:1 with the second most
massive luminous subhalo
| in the simulation volume (i.e., $\sim 2$ Mpc or up to $\sim
6~{}R_{\mathrm{vir}}$ for Latte, and approximately double for ELVIS).
Main branch | Star particles formed in the main halo (as defined by the simulation) before $t_{\mathrm{MR_{3:1}}}$
Building block | Star particles formed in luminous halos that join with the main branch halo before $t_{\mathrm{MR_{3:1}}}$
Proto-Milky Way | Star particles formed in the main branch and building block systems, that coalesce into a single galaxy before $t_{\mathrm{MR_{3:1}}}$
Figure 1: Left: Diagram of the merger tree of the m12b simulation in Latte up
to a lookback time of 8 Gyr. In this work, a proto-Milky Way is defined as the
amalgamation of the main branch halo (red) in a simulation plus all the
building blocks (other colours) that coalesce onto it before
$t_{\mathrm{MR_{3:1}}}$ (dashed magenta line, see Table 1 for details). Right:
Distance at which a star particle in the simulation is formed w.r.t. the
centre of the main host as a function of lookback time (i.e., age) for the
m12b simulation in Latte. Highlighted in red are the star particles associated
with the main branch, and as other colours the star particles associated with
the resolvable building blocks (i.e., luminous subhaloes) that join with the
main branch before $t_{\mathrm{MR_{3:1}}}$. By tracking each system over time,
we are able to identify the star particles associated with all the proto-Milky
Way fragments.

Before performing any analyses, it is imperative that we make clear our
definition of a proto-Milky Way. To do so, we set out to answer the following
two questions:

1) When does a proto-Milky Way form?

2) What constitutes a proto-Milky Way?

For this work, we define the time in which a system becomes a proto-Milky Way
when the main halo in the zoom-in region of the simulation box ($\sim 2$ Mpc
side box or up to $\sim$6 $R_{\mathrm{vir}}$ for Latte, and approximately
double for ELVIS) reaches a stellar mass ratio of 3:1 with the second most
luminous halo666This second most luminous halo is by definition not in our
sample of systems that comprise the proto-Milky Way.
($t_{\mathrm{MR_{3:1}}}$). We note that for the ELVIS galaxies, as there are
two main haloes in each simulation, we take the ratio of stellar mass between
the two main host haloes and the third most luminous subhalo within the
simulation. Albeit arbitrary, we argue that this definition is intuitively
sensible as: $i$) it defines a point in time in which a system is dominant in
its environment; $ii$) is quantifiable and repeatable across different
simulations. It has also been reasoned to be a soundproof way of pinpointing
the time in which a galaxy becomes dominant (Santistevan et al., 2020). 777We
note that Santistevan et al. (2020) also use another sensible choice, defined
as the time in which the main halo begins forming the majority of its present
day stars, signifying a dominant star formation episode of sustainable growth
from an in situ channel. Given our definition, it is then straight forward to
answer question 2); the proto-Milky Way is defined as the stellar population
resulting from the amalgamation between the main branch (i.e., the halo
tracing the formation of the main host in the simulation) and building blocks
(i.e., all the haloes that join with the main branch) before
$t_{\mathrm{MR_{3:1}}}$. Fig 1 shows a schematic of our definitions, where the
main branch is shown in red and building blocks as other colours. Moreover, a
list of definitions is provided in Table1; Table 2 displays the
$t_{\mathrm{MR_{3:1}}}$ times, the number of systems that coalesce to form a
proto-Milky Way (given our definitions), and the stellar/virial masses for
each proto-Milky Way system.

Table 2: List of approximate times in the simulation in which we define the proto-Milky Way to form ($t_{\mathrm{MR_{3:1}}}$), as well as the number of events (i.e., main branch+building blocks) that constitute the proto-Milky Way ($n_{\mathrm{systems}}$), and the stellar and virial mass of the proto-Milky Way at $t_{\mathrm{MR_{3:1}}}$. We note here that the number of building blocks we find is bounded by our choice to track only systems with 150 star particles or more (see Santistevan et al. (2020) for a comparison when tracking lower mass systems). We also list the corresponding references for each of the Latte/ELVIS simulated haloes: A. Hopkins et al. (2018b); B. Garrison-Kimmel et al. (2019a); C. Garrison-Kimmel et al. (2019b);, D. Garrison-Kimmel et al. (2017); E. Wetzel et al. (2016); F. Samuel et al. (2020). These (lookback) times should be interpreted as ages, where 13.8 Gyr is old and 0 Gyr is young. The average values for $n_{\mathrm{systems}}$ have been rounded to the closest integer. The $t_{\mathrm{MR_{3:1}}}$ agree well with the results found in Santistevan et al. (2020) (see their Table 1). Host |  $t_{\mathrm{MR_{3:1}}}$ [Gyr] |  $t_{\mathrm{MR_{3:1}}}$ [$z$] | $n_{\mathrm{systems}}$ | M${}_{\star_{3:1}}$ [$\times$109 M⊙] | M${}_{\mathrm{vir}_{3:1}}$ [$\times$1011 M⊙]
---|---|---|---|---|---
(isolated) |  |  |  |  |
m12bC | 9.76 | 1.59 | 8 | 8.80 | 4.43
m12cC | 9.05 | 1.32 | 5 | 3.56 | 4.93
m12fD | 12.15 | 3.77 | 4 | 0.38 | 2.86
m12iE | 11.80 | 3.18 | 6 | 0.35 | 1.12
m12mA | 9.94 | 1.68 | 5 | 3.69 | 5.00
m12rF | 8.05 | 1.02 | 7 | 4.29 | 2.61
m12wF | 11.60 | 2.92 | 9 | 0.35 | 1.19
(pairs) |  |  |  |  |
JulietC | 12.55 | 4.72 | 4 | 0.31 | 0.93
LouiseC | 11.90 | 3.32 | 4 | 0.99 | 1.02
RemusB | 7.90 | 0.98 | 9 | 20.14 | 6.41
RomeoC | 12.81 | 5.70 | 4 | 0.12 | 0.33
RomulusB | 10.00 | 1.70 | 6 | 6.85 | 3.40
ThelmaC | 12.40 | 4.29 | 2 | 0.11 | 0.78
Average (isolated) | 10.33 | 1.87 | 6 | 3.06 | 3.16
Average (pairs) | 11.26 | 2.56 | 5 | 4.75 | 2.14
Average (total) | 10.76 | 2.15 | 6 | 3.90 | 2.69

Furthermore, in order to track halos with a high enough number of star
particles to be resolvable in Latte/ELVIS, we choose to only track systems
with 150 star particles or more, using a similar technique to previous studies
(e.g., Necib et al., 2019; Panithanpaisal et al., 2021; Horta et al., 2023b).
More specifically, we track the evolution of each system (i.e., the star
particles) at every snapshot in the simulation over time until
$t_{\mathrm{MR}_{3:1}}$ or until the system no longer exists (i.e., it has
merged with the main branch system) using the halo catalogues, instead of
resorting to the merger trees. This allows us to identify all the star
particles in every halo, and to assign star particles to individual building
blocks. Our choice to only track systems with 150 star particles or more leads
to a minimum stellar mass of M⋆$\sim$1$\times$106M⊙ for subhaloes in Latte,
and M⋆$\sim$5$\times$105M⊙ for ELVIS, given the slight differences in particle
resolution between these two sets of simulations. For more information on the
limitations of this choice, see Section 6.

Lastly, in Fig 2 we show the stellar mass ratio between every main branch
(diamonds) and building block (circles) with the overall proto-Milky Way
population as a function of the stellar mass of each main branch/building
block. In this diagram, it is clear that there are two types of proto-Milky
Way populations: 1) proto-Milky Way systems formed by one clearly dominant
halo (namely, the main branch system. For example m12i). This is the case for
$\sim$60$\%$ (8/13) of the Milky Way-mass haloes in our sample; 2) proto-Milky
Way populations formed by two main systems (namely, the main branch and a
massive building block —$\gtrsim$1:5 mass ratio, see Fig 2—. For example,
m12m). This is the case for $\sim$40$\%$ of our sample (5/13 galaxies). In the
following sections we set out to study the chemical-kinematic properties of
these systems at $z=0$ to see if the different systems that form proto-Milky
Way populations can be distinguished.

Figure 2: Stellar mass ratio between each main branch (diamonds) or building
block (circles) and the proto-Milky Way (i.e., main branch$+$building blocks)
at $t_{\mathrm{MR_{3:1}}}$ as a function of the stellar mass of each main
branch/building block. There are five clear cases (namely, m12c, m12f, m12m,
m12w, and Romulus) which have a building block on the order of $\gtrsim$1:5
mass ratio with the proto-Milky Way population (i.e., blue shaded region).
This indicates that $\sim$40$\%$ of the thirteen proto-Milky Way systems
studied have a massive/dominant building block in addition to the dominant
main branch halo.

##  4 Results

In this section we consider each main branch and building block population
individually in order to examine how the mass, age, spatial distribution,
morphology, kinematics, and abundance patterns of these systems differ. We
also aim to examine the role of environment by comparing isolated Milky Way-
mass haloes with ones in Local Group environments.

###  4.1 Stellar mass and age

We begin by examining the stellar mass and the minimum age of a star particle
in each main branch and building block separately, shown in Fig 3. We choose
to study the minimum age instead of the mean age as we believe that it is a
better constraint for the time in which a system quenches star formation, and
thus ceases to evolve. The questions we aim to answer are: How old are the
systems that form a proto-Milky Way? How massive are the main branch and
building block systems that constitute a proto-Milky Way?

The left panel of Fig 3 shows the relationship between the minimum star
particle age and the stellar mass at $z=0$ of the main branch (diamonds) and
building blocks (circles) from each simulation. The middle and right panels
show histograms for stellar mass and minimum star particle age, respectively.
Here, we show the frequency of all systems together in black, isolated haloes
(i.e., Latte) in green, systems in Local Group environments (i.e., ELVIS) in
yellow, main branch systems in red, and building blocks in blue. We recall
that the resolution of the Latte and ELVIS suites are different (see Section 2
and the figure caption for details).

Figure 3: Left: minimum star particle age (i.e., the youngest star particle) in each main branch (diamond) and building block (circle) event as a function of its stellar mass for all thirteen Milky Way-mass haloes. Middle: histogram of the stellar mass of each system, split into environment and main branch/building block populations. The vertical dotted lines demark the minimum possible stellar mass of a halo in the simulation given our choice to track luminous subhalos with 150 star particles or more and the simulations resolution, for Latte (red) and and ELVIS (blue). Right: histogram of the minimum star particle age, illustrated as in the middle panel. Table 3: Summary of the mean, median, and 1-$\sigma$ values for the stellar mass and minimum star particle age for our sample of events comprising proto-Milky Way systems. sample |  $\langle$log10(M${}_{\star})$$\rangle$ [M⊙] |  med(log10(M⋆)) [M⊙] |  $\sigma_{\log_{10}(\mathrm{M}_{*})}$ [M⊙] |  $\langle$agemin$\rangle$ [Gyr] |  med(agemin) [Gyr] |  $\sigma_{\mathrm{age_{min}}}$ [Gyr]
---|---|---|---|---|---|---
proto-Milky Way (all) | 7.39 | 7.08 | 1.17 | 11.36 | 11.64 | 1.16
Environmental difference |  |  |  |  |  |
proto-galaxy (pairs) | 7.21 | 6.79 | 1.25 | 11.50 | 11.97 | 1.41
proto-galaxy (isolated) | 7.50 | 7.31 | 1.09 | 11.26 | 11.51 | 0.94
Classification difference |  |  |  |  |  |
main branch | 9.09 | 8.99 | 0.71 | 10.76 | 11.60 | 1.66
building block | 7.02 | 6.81 | 0.88 | 11.10 | 11.20 | 1.31
All |  |  |  |  |  |
main branch (pairs) | 8.96 | 8.74 | 0.85 | 11.26 | 12.15 | 1.76
main branch (isolated) | 9.19 | 9.55 | 0.56 | 10.33 | 9.94 | 1.43
building block (pairs) | 6.75 | 6.64 | 0.88 | 11.57 | 11.97 | 1.29
building block (isolated) | 7.178 | 7.05 | 0.84 | 11.43 | 11.52 | 0.68

From inspection of the middle and right panels of Fig 3, we find that, overall
(black solid line), the mean stellar mass of the subhaloes that comprise
proto-Milky Way systems is $\langle$M⋆$\rangle$$\sim$4$\times$107M⊙, and the
mean minimum age is $\langle$agemin$\rangle$$\sim$11.4 Gyr (or
$\langle$$z$$\rangle$$\sim$2.8). However, the range in these parameters is
significantly large, and for both cases, we find that the distribution is
skewed (or for the case of the stellar mass, even possibly bimodal). Given our
set of choices and assumptions, we find that the stellar masses for all these
systems ranges from $\sim$5$\times$105$<$M⋆$<$2$\times$1010M⊙, and the minimum
age can reach anything between $\sim$8$<$agemin$<$13 Gyr (or
$\sim$0.7$<$$z$$<$6.5). Furthermore, when splitting the distributions by
environment and comparing the haloes that live in isolated environments from
those that live in a Local Group-like environment, we find that the mean
values of the stellar masses of systems in Local Group-like environments is
slightly lower than for systems in isolation. The difference in the median
stellar mass is on the order of half a magnitude. Furthermore, there are also
some slight differences in the distribution of the minimum star particle age,
whereby the pairs tend to favour older populations. Statistically, the mean
values show no clear differences ($\langle$agemin,pairs$\rangle$ = 11.5 Gyr
and $\langle$agemin,isolated$\rangle$ = 11.26 Gyr). However, the median
differs by approximately 0.5 Gyr for both the stellar mass and minimum age
(see Table 3).

When comparing the main branch populations (red) with their building block
counterparts (blue), we see substantial differences in their stellar masses.
Specifically, we find that main branch systems tend to be more massive,
approximately two orders of magnitude larger (i.e., $\langle$M⋆,mb$\rangle$ =
1$\times$109M⊙ and $\langle$M⋆,bb$\rangle$ = 1$\times$107M⊙). However, we note
that some building blocks do reach high stellar masses, on the order of
M⋆$\sim$2$\times$109M⊙, making the range in the distribution for building
blocks much larger than for main branch debris. We see similar differences in
the mean and spread of the minimum star particle age between main branch and
building blocks debris as we did when comparing environment, on the order of
$\sim$0.4 Gyr.

In the left panel of Fig 3 we show the distribution of minimum star particle
age (i.e., the youngest) as a function of stellar mass for each event in each
Milky Way-mass halo separately, where here we distinguish main branch
populations (diamonds) and building blocks (circles). As can be seen from this
illustration, there is a relation between the stellar mass of a system and the
minimum star particle age, whereby a system’s minimum stellar age decreases
with increasing stellar mass. This is likely because those systems that are
forming stars for longer are able to build up more mass (a similar argument as
for the mass-metallicity relation of galaxies, see also Figure 3 of Cunningham
et al., 2022). These results from this section are summarised in Table 3.

In summary, when examining the stellar masses and minimum star particle ages
(a proxy for quenching/accretion time) of the events that comprise proto-
galaxy systems, we have found the following:

  * •

Overall, proto-Milky Way populations are largely formed from a range of
different mass and age systems, but are predominantly comprised by one or two
major systems of similar mass to the LMC (i.e., M⋆$\sim$1$\times$109M⊙) and
3–5 smaller mass building blocks (i.e., M⋆$\sim$4$\times$107M⊙), see Fig 2 and
Fig 3 (see also Fig 4 in Santistevan et al., 2020). The massive systems
typically arise from the main branch progenitor, and in the case where there
are two dominant systems, also from a massive building block. We note that the
number of building blocks we find are subject to our choice to only study
systems of a minimum stellar mass of M⋆$\sim$1$\times$106M⊙ for Latte, and
M⋆$\sim$5$\times$105M⊙ for ELVIS. Santistevan et al. (2020) showed that when
lowering this limit to lower masses ($\sim$105 M⊙) and looking at all subhalos
that coalesce with the main branch before $z=0$, this number grows
significantly to $\sim$100 building blocks.

  * •

There is a small difference when comparing the mean values of haloes in pairs
and in isolation. Here, pairs tend to be less massive (med(M⋆) $\sim$
8$\times$106) and older (med(agemin) $\sim$ 12 Gyr) than systems in isolation
(med(M⋆) $\sim$ 3$\times$107 and med(agemin) $\sim$ 11.5 Gyr). This is likely
because systems in pairs assemble earlier (Santistevan et al., 2020).

  * •

For this suite of simulations and given our definitions, we find a clear
difference in the mean stellar mass value of main branch populations when
compared to building blocks. More specifically, main branch populations are
typically two orders of magnitude more massive
($\langle$M⋆,mb$\rangle$$\sim$1$\times$109M⊙) when compared to building blocks
($\langle$M⋆,bb$\rangle$$\sim$1$\times$107M⊙). Conversely, this difference is
less pronounced in the mean values of the minimum star particle age when
comparing main branch debris to their counterpart building blocks, showing
only a difference of $\sim$0.4 Gyr.

###  4.2 Spatial distribution and morphology

Figure 4: Total stellar mass fraction as a function of spherical radius (at
$z=0$) for all the main branch (dashed) and building block (solid) events in
each Milky Way-mass halo. Each profile is colour coded by the systems’
respective stellar mass. The vertical dashed lines mark 5 kpc (cyan) and 20
kpc (gray). The majority of the proto-Milky Way material (i.e., main branch
and massive building blocks) are contained within a small spatial volume,
close to the host’s centre. However, there is a range of spatial profiles for
lower mass building blocks.

We now set out to examine the spatial distribution of the different
contributors to a proto-Milky Way. In doing so, we aim to tackle the following
questions: Where is most of the mass contained in the main host at $z=0$? What
is the morphology of the debris at $z=0$?

Fig 4 characterizes where mass is deposited by showing the total (cumulative)
mass fraction of the debris of every main branch (dashed) and building block
(solid) event comprising the proto-Milky Way systems as a function of present
day spherical radius from the centre of the Milky Way-mass host. Each main
branch/building block event is colour coded by their respective stellar mass.
To guide the readers eye, we also plot a vertical line at 5 kpc (cyan) and 20
kpc (grey).

Figure 5: Histograms of the 50th (top) and 95th (bottom) percentiles of the
spherical radii for all (black), isolated haloes (green), haloes in pairs
(yellow), main branch systems (red), and building blocks (blue). The majority
of the stellar populations of proto-Milky Ways are contained within $\sim$40
kpc, with 50$\%$ of the mass contained within $\sim$5-10 kpc from the host’s
centre at present day. Figure 6: Ellipticity (i.e., the ratio of 95th
percentiles between cylindrical height and radius —i.e., the aspect ratio—
subtracted from unity) as a function of the stellar mass for main branch
debris (diamonds) and building block events (circles). More massive systems
tend to adopt a more oblate shape.

As is to be expected, there is a range of spatial distributions for all the
main branch and building block events. However, the majority of the mass from
these systems is predominantly contained within $\sim$10 kpc from the centre
of the host halo. This is especially the case when examining either main
branch progenitors (dashed lines) and/or more massive building block debris
(darker solid lines). For all Milky Way-mass galaxies the main branch system
contains $\sim$50$\%$ of its mass within $\sim$10 kpc from the host’s centre,
and $\sim$95$\%$ within $\sim 30-40$ kpc. For the case of the building blocks,
the most massive events (M⋆$\gtrsim$5$\times$108 M⊙) typically contain
$\sim$50$\%$ of their mass within $\sim 5-20$ kpc, and $\sim$95$\%$ of their
mass within $\sim 20-40$ kpc, in a similar fashion to the main branch
population. This is especially the case for m12f, m12m, and Romulus, where the
main branch progenitor follows a very similar spatial profile to the most
massive building block in that system, making these almost indistinguishable
spatially. For the lower mass building blocks (M⋆$\lesssim$5$\times$106 M⊙),
we see a wide range of profiles, ranging from $\sim 10-50$ kpc for 50$\%$ of
their enclosed mass and $\sim 20-100$ kpc for 95$\%$ of their enclosed mass.
We find a small difference in the spatial distribution of the debris of an
event that occurs in a halo that is isolated versus a halo that is in a Local
Group environment, whereby the latter tend to present proto-Milky Way
populations more concentrated toward the host’s centre. These results can be
more easily digestible in Fig 5, where we show histograms of $r_{50}$ and
$r_{95}$ (i.e., the 50th and 95th percentiles of the spherical radii for each
system), following the colour coding from Fig 3.

Fig 6 provides a summary of morphology by showing the aspect ratio ($c/a$)
subtracted from unity, as a function of their stellar mass. As in previous
figures, the main branch populations are shown as diamonds and building block
systems as circles. Given that the 95th percentiles represent almost the
extent of the distribution, we use these values to define two axes: a semi-
major one ($a=R_{95}$) and a semi-minor one ($c=z_{95}$), both in units of
kpc. We define these quantities so that we can measure the aspect ratio
($c/a$) and in turn their ellipticity, defined as $\varepsilon$ = 1 – $c/a$.
Given this definition, a value of $\varepsilon$ = 1 demarks a perfect circle,
and an $\varepsilon$ value closer to 0 corresponds to a more squashed ellipse.
We note that these aspect ratios are defined with respect to the principal
axis of the host at $z=$0.

We find that the there is a dependence on $\varepsilon$ with stellar mass,
whereby more massive systems tend to adopt a more oblate distribution at
present day when compared to their lower mass counterparts, which adopt a
wider range of morphologies. As we saw in Section 4.1, the main branch debris
are typically the higher mass events. This leads to conclude that main branch
progenitors adopt a more oblate morphology when compared to the lower mass
building block counterparts. However, we find that the more massive building
block debris also adopt an oblate distribution, following the main branch
progenitors.

In summary, our investigation into the spatial distribution and morphology of
proto-Milky Way populations and early building blocks have led to the
following results:

  * •

For all host haloes the main branch debris contain 50$\%$ of their mass within
$\sim$10 kpc from the host’s centre, and 95$\%$ within $\sim 40-60$ kpc (Fig 4
and Fig 5). This is also the case for the most massive building blocks
(M⋆$>$5$\times$108 M⊙). Conversely, for lower mass building blocks, there is a
wide range of spatial distributions. The value of 50$\%$ of the total mass can
range from $\sim 10$ kpc to $\sim 100$ kpc. While there are differences in the
spatial distribution of main branch debris and building block counterparts, at
fixed mass distinguishing these is challenging.

  * •

There is a small difference in the spatial distribution of debris in isolated
hosts and haloes in a Local Group-like environments, whereby pairs tend to
present proto-Milky Way populations more concentrated towards the host’s
centre (Fig 5).

  * •

There is a relation between an event’s stellar mass and present day
morphology, whereby more massive systems tend to be more oblate in shape when
compared to lower mass counterparts. This is true for both main branch
progenitors and massive building blocks. At lower mass, there is a larger
spread in $\varepsilon$ values, indicating a wider range of possible present
day morphologies, including more spherical distributions.

###  4.3 Kinematic properties

In this section, we set out to examine the kinematic properties of populations
comprising proto-Milky Way systems. Specifically, we aim to answer the
following questions: How rotationally supported are the main branch and
building block systems that form the proto-Milky Way at present day? Are
proto-Milky Way populations rotating on prograde or retrograde orbits?

To quantify the amount of rotational support in each system, we compute the
ratio of the rotational velocity over the total velocity for all star
particles in a given population,
$\kappa_{\mathrm{rot}}$/$\kappa_{\mathrm{tot}}$. These two velocity quantities
are defined in the coordinate system centred on the host Milky Way-mass galaxy
at $z=0$, and are computed by finding the centre of the Milky Way-mass host
using its baryonic particles at that redshift. This can be estimated by taking
the ratio of the kinetic energy in the rotational direction over the total
kinetic energy (McCarthy et al., 2012), such that:

$\kappa_{\mathrm{rot}}=(L_{z}/R)^{2}/2,$ (1)

assuming $R$ = $\sqrt{X^{2}+Y^{2}}$, $L_{z}=|\vec{R}\times\vec{v_{\phi}}|$,
and

$\kappa_{\mathrm{tot}}=v_{\mathrm{tot}}^{2}/2,$ (2)

where

$v_{\mathrm{tot}}=\sqrt{v_{X}^{2}+v_{Y}^{2}+v_{Z}^{2}}.$ (3)

We normalise this quantity by taking the ratio of this value with the amount
of rotational support determined in the young (age$<4$ Gyr) disc
($d_{\mathrm{form}}<30$ kpc) of the host of each Milky Way-mass halo,
$\kappa/\kappa_{\mathrm{disc}}$ (where
$\kappa=\kappa_{\mathrm{rot}}$/$\kappa_{\mathrm{tot}}$).

Fig 7 shows the mean value of $\kappa$/$\kappa_{\mathrm{disc}}$ for all the
main branch (diamonds) and building block (circles) systems as a function of
their mean ellipticity. Milky Way-mass haloes in isolation (pairs) are shown
as full (empty) markers. A value of $\kappa$/$\kappa_{\mathrm{disc}}=1$
signifies that a system is as rotationally supported as the young disc in that
simulation. A value of $\kappa$/$\kappa_{\mathrm{disc}}>1$ implies that it is
more rotationally supported, and a value of
$\kappa$/$\kappa_{\mathrm{disc}}<1$ that it is less rotationally supported. In
right panel, we plot the mean azimuthal velocity of every proto-Milky Way
population.

We find that the majority of the debris from both the main branch and building
block systems have a mean $\kappa$/$\kappa_{\mathrm{disc}}$ value between 0.25
and 0.8, indicating not a strong amount of rotational support in the system
when compared to the host’s young disc. However, their rotational support is
not 0, indicating that these systems are not in a purely isotropic
distribution (Santistevan et al., 2022; Yu et al., 2023). In fact, their mean
azimuthal velocities show that the majority of these systems have
predominantly prograde orbits ($\langle v_{\phi}$$\rangle>$0 km/s). This means
that they are rotating in the same direction as the disc, although at a much
slower pace. Our results suggest that all these systems, either defined as
main branch progenitor or building block, show weak but systematic net
rotation in the plane of the host’s disc at $z=0$ despite many of them
appearing flattened/oblate in morphology (see Section 4.2).

To summarise, our results for the kinematic profiles reveal that:

  * •

Despite the range of morphologies, the majority of proto-Milky Way populations
show weak but systematic net rotation in the plane of the host’s disc at $z=0$
(i.e., 0.25$<$$\kappa$/$\kappa_{\mathrm{disc}}$$<$0.8). This indicates that
proto-Milky Way systems are not in a purely isotropic distribution.

  * •

The majority of proto-Milky Way populations have prograde orbits ($\langle
v_{\phi}$$\rangle>$0).

  * •

We find no dependence of the mean $\kappa$/$\kappa_{\mathrm{disc}}$ value of
proto-Milky Way populations with environment (see the left panel of Fig 7).

Figure 7: Left: mean of the rotational support of each proto-Milky Way
population normalised by the rotational support of its host’s disc
($\kappa/\kappa_{\mathrm{disc}}$) as a function of ellipticity. Rotational
support is defined as the ratio of rotational energy ($\kappa_{\mathrm{rot}}$)
in the direction of the host’s galactic disc at present day and the total
kinetic energy ($\kappa_{\mathrm{tot}}$). The disc population is comprised by
star particles with age$<$4 Gyr formed in the main host,
$d_{\mathrm{form}}<$30 kpc. Right: mean azimuthal velocity in the direction of
the host’s disc at present day as a function of ellipticity. Systems in
isolation (pairs) are shown as full (empty) markers. The majority of main
branch and building block debris show low level of systematic net rotation
(0.25$\lesssim$$\langle$$\kappa$/$\kappa_{\mathrm{disc}}$$\rangle$$\lesssim$0.8)
when compared to their host present day disc. However, this rotational support
is not zero. In fact, the majority of proto-Milky Way populations show
prograde motion (0$\lesssim$$\langle v_{\phi}$$\rangle\lesssim 50$ km/s), that
can reach up to $\langle v_{\phi}$$\rangle\sim 100-150$ km/s (m12i).

###  4.4 Chemical compositions

In this Section, we examine what is likely to be the most pristine observable
tracer of stellar populations: their chemical compositions. While FIRE tracks
eleven different elemental species, they trace three nucleosynthetic channels:
contributions from core-collapse supernovae, (CCSN), SN type Ia, and stellar
winds. We set out to explore the distribution of the debris of these events in
the chemical plane tracing the contribution of CCSN and SNIa, choosing Mg as
our $\alpha$ tracer. The main question we aim to tackle in this Section is:
what differences in the chemical compositions (and thus, the star formation
histories) can we expect from main branch progenitors when compared to their
building block counterparts?

Fig 8 shows the median [Mg/Fe] and [Fe/H] values for building blocks (circles)
and main branch debris (diamonds), where each system is colour coded by their
mean stellar mass (left) and minimum star particle age (right). Empty markers
correspond to systems in Local Group-like environments (i.e., the ELVIS suite)
and filled markers correspond to events in isolated haloes (i.e., the Latte
suite).

Fig 8 demonstrates that, on average, the main branches are more metal-rich
than the building blocks. This is a natural consequence of the mass-
metallicity relation (Kirby et al., 2013, 2020, e.g.,). This fact could be
leveraged to disentangle the lower mass building blocks from their main
branches using their abundances (e.g., employing methods presented in
Cunningham et al. 2022, Horta et al. 2023a, Deason et al. 2023, for example).
However, some of the more massive building blocks have chemical compositions
consistent with the main branches. As a result, it could be challenging to
disentangle populations from the main branch from populations from the most
massive building block when a proto-Milky Way has two significant contributors
(which is the case for roughly $\sim$40$\%$ of our Milky Way-mass galaxy
sample). We investigate this further in Fig 9, in which we compare the
metallicity distribution functions (MDFs) for the main branch (dashed) and
most massive building block (solid) in both m12f (black) and m12m (red). In
the case of m12f, the mean of the main branch MDF is $\sim$0.1–0.2 dex higher
compared to the most massive building block. Conversely, the main progenitor
in m12m is identical in [Fe/H] when compared to its massive building block
counterpart. These results show that, quantitatively, it would be extremely
difficult to distinguish massive building blocks from their main branch system
counterparts using their MDFs.

When comparing the debris from host haloes that evolve in different
environments, we find that there are possible subtle differences.
Specifically, we find that the main branch and building block events in host
haloes simulated in Local Group-like environments (pairs) present on average
slightly higher median [Mg/Fe] and/or [Fe/H] values when compared to isolated
hosts. This result is consistent with our finding that systems in Local Group-
like environments evolve either faster and/or earlier, when compared to
systems in isolated hosts (Santistevan et al., 2020). However, within the
spread of the distributions (i.e., the uncertainties in Fig 8), which is on
average $\sim 0.05-0.1$ dex for [Mg/Fe] and $\sim 0.3$ dex for [Fe/H], we note
that this result is tentative.

Figure 8: Median [Mg/Fe] and [Fe/H] values for main branch (diamonds) and
building block (circles) populations for isolated Milky Way-mass haloes
(filled markers) and those in pairs (empty markers), colour coded by stellar
mass (left) and minimum star particle age (right). The average 1-$\sigma$ is
0.33 dex for [Fe/H], and 0.08 dex for [Mg/Fe]. On average, main branch
progenitors and the most massive building blocks present more enriched
chemical compositions when compared to lower mass building blocks; this is a
result of the galaxy mass-metallicity relation. Figure 9: Histograms of the
metallicity distribution functions (MDF) of main branch progenitors (dashed)
and the most massive building blocks (solid) events in the m12f (black) and
m12m (red) simulations. The MDF profiles are qualitatively the same for these
populations. This result highlights the difficulty of distinguishing main
branch from massive building block populations based on their MDFs.

In summary, our findings on the chemical properties of main branch and
building block events can be summarised as follows:

  * •

On average, the [Mg/Fe]-[Fe/H] distribution of main branch and building block
events shows a strong dependence on the debris mass and age. On average, the
main branch systems tend to show more enriched [Mg/Fe]-[Fe/H] chemistry than
the majority of their building block counterparts. Here, main branch
progenitors either have a higher [Mg/Fe] value at fixed [Fe/H], or higher
[Fe/H] overall. However, we find that this is also the case for the most
massive (and youngest) building blocks. As a result, it is difficult to
distinguish, qualitatively, these two types of populations based on the
average [$\alpha$/Fe]-[Fe/H] abundances alone. More detailed chemical
abundance information from elemental species synthesised in more exotic
nucleosynthetic channels, which we have not been able to examine in this work
(but see Horta et al., 2023a), may hold the clues to disentangling these
dominant proto-Milky Way fragments.

##  5 Discussion

###  5.1 What constitutes a proto-Milky Way?

In this work, we have defined a proto-Milky Way to be comprised by a main
branch and building blocks systems (Table 1 and Fig 1). Given our definitions,
our findings suggest that on average the proto-galaxy of a Milky Way-mass halo
is comprised of five to six systems with M${}_{\star}\gtrsim 5^{5}$M⊙.
888However, these results are subject to our choice to only study systems with
150 star particles. If one was to lower this mass limit and was to look at all
systems that coalesce with the main branch before $z=0$ (e.g., Santistevan et
al., 2020), this number grows significantly, up to $\sim$100 building blocks..
Although likely insignificant, it is worth noticing that the average number of
systems that form a proto-Milky Way is smaller in the Local Group Milky Way-
mass galaxy simulations (5) when compared to the isolated haloes (6). This is
likely because, as found in Fig 3, Milky Way-mass galaxies in Local Group
environments assemble earlier than those in isolation (Santistevan et al.,
2020).

Given this finding, it is interesting to ask what fraction of the proto-
galaxy’s (stellar) mass is contributed by the main branch/building block
progenitors? To answer this question, in Fig 10 we show the mass difference
between each main branch progenitor and its counterpart building blocks as a
function of the mass of the main branch progenitor. We find that the majority
of the building blocks are of lower mass when compared to the main branch
progenitor, on the order of one to four orders of magnitude difference, with a
stellar mass ranging between $\gtrsim$ 5$\times$105$<$M⋆$<$1$\times$108 M⊙
(see also Fig 2). This is clearly the case for eight out of the thirteen Milky
Way-mass haloes (namely, m12b, m12i, m12r, Remus, Thelma, Louise, Romeo, and
Juliet), and would indicate that the majority of these proto-Milky Way systems
are comprised by one dominant progenitor and several ($\sim 4-5$) smaller mass
building blocks. Conversely, for the other five proto-Milky Way systems, we
find that there is at least one building block that has a mass ratio with the
main branch progenitor that is within one order of magnitude (i.e., has a
stellar mass greater than $\sim$5$\times$108 M⊙). For these haloes (namely
m12c, m12f, m12m, m12w, and Romulus), the proto-Milky Way is formed by two
dominant systems, and a smaller number ($\sim 3-4$) of lower mass building
blocks. This result is interesting, as it suggests that proto-galaxy systems
in Milky Way-mass galaxies can be formed by either one or two main subhalos at
high redshift, in addition to a small number of lower mass building blocks. We
find that the mass ratios between main progenitor branches and building blocks
does not show a correlation with either the main branch progenitor’s stellar
mass or environment.

###  5.2 When does a proto-Milky Way form?

Table 2 shows that there is a range of times in which a proto-Milky Way forms,
given our definition of when a proto-Milky Way emerges (see Section 3). These
values range from as late as $t_{\mathrm{MR_{3:1}}}$ = 8.05 Gyr ($z$=1.02) to
as early as $t_{\mathrm{MR_{3:1}}}$ = 12.9 Gyr ($z$=6.08). We find that
overall the proto-Milky Way systems can be grouped into three main camps: an
early forming group, an intermediate forming group, a late forming one. In the
early forming group we see m12f, m12i, m12w, Thelma, Louise, Romeo, and
Juliet, which all reveal a value of $t_{\mathrm{MR_{3:1}}}$ $\gtrsim$ 11.5 Gyr
($z$ $\gtrsim$ 2.8). In the intermediate age forming group we see m12m, and
Romulus, that have 10 $\lesssim$ $t_{\mathrm{MR_{3:1}}}$ $\lesssim$ 11.5 Gyr
(1.7 $\lesssim$ $z$ $\lesssim$ 2.8). In the late forming group, we see m12b,
m12c, m12r, and Remus, which all show 8 $\lesssim$ $t_{\mathrm{MR_{3:1}}}$
$\lesssim$ 10 Gyr ($z$ $\lesssim$ 1.7). It is interesting to note that the
members of the intermediate group (namely, m12m, and Romulus) are the majority
of the haloes that have the smallest mass fraction between the main progenitor
branch and the most massive building block (see Fig 10). The difference in
formation times of these Milky Way-mass galaxies has also been divided into
three similar groups based on the time in which their stellar discs settle
(McCluskey et al., 2023), and based on the transition from spheroids to thick
and thin discs (Yu et al., 2023).

These results are consistent with the results from Santistevan et al. (2020))
when they use the same definition as used here ($t_{\mathrm{MR_{3:1}}}$); for
the emergence of the Milky Way (though we note that they find that Milky Way-
mass galaxies emerge around $z\sim$3–4 when using a definition based on in-
situ star formation).

Figure 10: Ratio of the stellar mass between main branch progenitors and their
building block counterparts as a function of the main branch progenitor’s
stellar mass for each Milky Way-mass simulation at $t_{\mathrm{MR}_{3:1}}$.
Five out of thirteen proto-Milky Way’s host a building block of similar mass
to the main branch progenitor (blue shaded region).

###  5.3 The role of environment

When comparing the properties of systems in isolation and in pairs, we have
found a few subtle differences. For example, the median stellar mass per event
is lower for systems in pairs ($\sim$8$\times$106 M⊙) when compared to
isolated ($\sim$3$\times$107 M⊙). Similarly, the minimum star particle age is
$\sim$12 Gyr for pairs, and $\sim$11.5 for isolated haloes. This implies that
systems in Local Group environments assemble more rapidly and from lower mass
events when compared to systems in isolation, as pointed out in previous work
(Santistevan et al., 2020).

In terms of the spatial, kinematic, and chemical properties, we have found
subtle differences with environment. For example, the radial distribution for
50$\%$ of the mass for systems in pairs is more centrally concentrated towards
the host’s centre than systems in isolation (see Fig 5). Another example is
that we found that systems in pairs tend to present more enriched
[$\alpha$/Fe] and/ or [Fe/H] values than systems in isolation (Fig 8).

###  5.4 Is it possible to distinguish the different systems that form a
proto-Milky Way from one another?

Figure 11: Left: Density level contours of the X-Z positions (i.e., edge-on
projections) at $z=0$ for star particles belonging to the main branch system
(red), most massive building block (black), the old disc (blue), and the young
disc (green) in m12i. Here, the old(young) disc is defined as star particles
formed in the main host halo between
4$<$$\tau$[Gyr]$<$$t_{\mathrm{MR_{3:1}}}$(4$<$$\tau$[Gyr]). Middle: Same as
left, but now in the cylindrical Toomre diagram. Right: The same as left, but
now in the [Mg/Fe]-[Fe/H] plane. m12i and m12m differ in the fact that the
former proto-Milky Way is comprised primarily by one main branch system
($\sim$90$\%$ of the stellar mass) and m12m is comprised by two systems, a
main branch system ($\sim$65$\%$ of M⋆) and a massive building block
($\sim$35$\%$ of M⋆). The kinematic distributions of main branch, building
block, and old disc are qualitatively similar and are different from the young
disc. This is more pronounced in m12i than m12m. The [Mg/Fe]-[Fe/H]
compositions of the main branch and building block system are also extremely
similar, but are different from the young disc. The main branch and building
block overlap with the metal-poor old-disc; this overlap is bigger for m12m
than or m12i.

The main aim of this work was to identify the different systems that comprise
proto-Milky Way populations and to compare their observable properties.
Specifically, we have examined their masses, ages, spatial distribution,
morphology, kinematics, and chemical abundances. In this section, we summarise
these findings and put them in context with the latest work.

In terms of the mass and age, in Section 4.1 we show that overall there are
clear differences between the average main branch progenitor and their
counterpart building block populations. Main branch systems tend to be two
orders of magnitude more massive in stellar mass
($\langle$M⋆,mb$\rangle$$\sim$1$\times$109M⊙ vs
$\langle$M⋆,bb$\rangle$$\sim$1$\times$107M⊙). However, there is no clear
difference between their respective minimum star particle age (i.e., quenching
time) distributions. Although the mean values of stellar mass differ between
main branch and building block systems overall, for five out of the thirteen
Milky Way-mass haloes ($\sim$40$\%$), the mass ratio between the main branch
progenitor and the most massive building block is below one order of magnitude
(see Section 5.1 and Fig 10). In these cases, the massive building block
population contributes between $\sim$1/3 to $\sim$1/5 of the proto-Milky Way
stellar mass (see Fig 2). 999We recall that the maximum mass ratio between a
building block and a main branch system is limited by our definition of when a
proto-Milky Way forms ($t_{\mathrm{MR}_{3:1}}$, see Table 1), so it is
possible that building blocks can reach an even higher mass ratios.. This
implies that for over a third of the Milky Way-mass galaxies, it is two
systems that dominate the mass of the resulting proto-Milky Way population.
This finding highlights a possible dichotomy in the formation of proto-Milky
Way systems. On the one hand, a proto-Milky Way can form from one dominant
main branch progenitor plus a handful of additional lower mass building
blocks. On the other, a proto-Milky Way can form from two main haloes (i.e.,
the main branch and a massive building block) and 3–4 lower mass building
blocks.

In Section 4.2 we set out to examine the radial profiles of main branch
progenitors and building blocks (see Fig 4 and Fig 5). Our results showed that
the majority of the mass (95$\%$) of main branch systems is contained within
40–60 kpc from the host’s centre at present day, and 50$\%$ of this mass is
within $\sim 5-10$ kpc. On the other hand, building blocks show wider variety,
ranging from 5–50 kpc for $\sim 50\%$ of their mass. However, we do note that
for the most massive building blocks (M⋆$\gtrsim$5$\times$108M⊙), the spatial
distribution is largely similar to the main branch progenitors, especially in
the case of m12f, m12m, and Romulus. This result suggests that the innermost
regions of the Milky Way-mass galaxies is where you would expect to find the
oldest populations in the Galaxy (El-Badry et al., 2018; Fragkoudi et al.,
2020) that comprise the proto-Milky Way, arising from both a main branch
progenitor and (possibly) the most massive building blocks (Horta et al.,
2021a; Rix et al., 2022). Due to the strong spatial overlap between main
branch systems and massive building blocks, it would be difficult to
distinguish these based on spatial distribution information alone. Along those
lines, we also found that it would be difficult to distinguish these systems
based on their morphology (Fig 6), as we found both main branch systems and
massive building blocks to have a similar oblate shape for a given stellar
mass. However, distinguishing the dominant fragments of proto-Milky Way
systems from their lower mass counterparts, based on their spatial
distribution and morphology profiles of full samples, in principle could be
possible.

We also set out to examine the kinematics of main branch progenitors and
building blocks in Section 4.3. Here, we attempted to see if there were any
signatures of rotation in this population (Fig 7), and to assess if there is a
preference for prograde or retrograde orbits. Our results show that both main
branch progenitors and building block systems show weak but systematic net
rotation in the plane of the host’s disc at $z=0$ (i.e.,
0.25$\lesssim$$\langle$$\kappa/\kappa_{\mathrm{disc}}$$\rangle$$\lesssim$0.8).
This indicates that their distribution is not purely isotropic. Moreover, we
also found that the majority of the populations rotate on prograde orbits
($v_{\phi}>0$ km/s). It is important to note that recent studies have reported
that Milky Way-mass galaxies in FIRE experience different phases of growth (Yu
et al., 2021; Gurvich et al., 2023; Hopkins et al., 2023; McCluskey et al.,
2023; Semenov et al., 2023; Yu et al., 2023), which can be characterised based
on differences in their stellar kinematics and amount of rotational support.
Yu et al. (2023) find that the orbits of star particles formed in the main
branch spheroid at early times are more radial than those formed at later
times in a more settled disc (see also Gurvich et al., 2023 and McCluskey et
al., 2023). This result is related to our finding, and suggests that
populations formed at early times are likely to not be on as rotationally
supported orbits when compared to populations formed later in settled discs.

Lastly, in Section 4.4 we set out to compare the chemical compositions. Fig 8
shows the median [Mg/Fe]-[Fe/H] values for main branch progenitors and
building blocks. Overall, we find a similar result as when examining the
spatial and kinematic properties. On average, main branch progenitors present
a more enriched [Mg/Fe]-[Fe/H] value, either in the form of more [Fe/H]
overall or a higher [Mg/Fe] ratio for fixed [Fe/H], when compared to the bulk
of building block populations. As argued in previous studies (Cunningham et
al., 2022; Horta et al., 2023b, e.g.,), this is caused by the effect of the
mass-metallicity relation of galaxies. More massive systems that are able to
form stars for longer enrich their interstellar medium with more metals and
evolve chemically faster than lower mass systems that quench star formation
earlier (Kirby et al., 2013). This is the reason why the bulk of the main
branch systems (that are more massive and evolve for longer) sit towards the
metal-rich end of Fig 8 when compared to lower mass building blocks. However,
as discussed at the end of Section 4.4, the more massive building blocks
resemble qualitatively main branch progenitors in the [Mg/Fe]-[Fe/H] plane. To
further analyse these specific cases, we took two case examples (m12f and
m12m) where the most massive building block is of similar mass to the main
branch progenitor (see Fig 10) and examined their metallicity distribution
functions (MDF; Fig 9). In doing so, we found a similar result. We find that
their MDF profiles are qualitatively, and quantitatively, similar. This result
stresses how difficult it is to distinguish, on the basis of [Fe/H]
information alone, main branch populations from its corresponding massive
building block counterparts. However, distinguishing the dominant fragments
from lower mass ones is likely possible (Deason et al., 2023).

In summary, our results reveal that it is likely possible to distinguish
smaller mass building blocks from the more massive main branch progenitor
systems in chemical-kinematic samples. However, when attempting to distinguish
main branch progenitors from the more massive building blocks, this task will
likely result rather complicated, as they appear indistinguishable in almost
every chemical and kinematic plane studied. Despite this challenging task, we
argue that obtaining detailed elemental abundance ratios for a wide spectrum
of elements (probing different nucleosynthetic channels) for large samples of
metal-poor stars in the inner Galactic regions of the Milky Way may result
profitable (e.g., Kollmeier et al., 2017), as it is likely that these stellar
tags are the only possible diagnostic to distinguish such populations.

###  5.5 Can we distinguish if a proto-Milky way formed from one or two
dominant systems?

Fig 11 shows the X-Z positions, Toomre diagram distribution, and
[Mg/Fe]-[Fe/H] density contour distributions for the main branch (red), the
most massive building block (black), the old disc (blue), and young disc
(green) in the m12i and m12m simulations. Disc stars are selected to be formed
within 30 kpc of the main host after $t_{\mathrm{MR_{3:1}}}$, where the old
(young) disc is older (younger) than 4 Gyr. We choose to compare m12i with
m12m, as m12i is a clear case where the proto-Milky Way formed from one
dominant main branch system (M${}_{\star}=$ 3.51$\times$108M⊙), and m12m is a
clear example of a proto-Milky Way formed from two dominant systems (a main
branch of mass M${}_{\star}=$ 3.69$\times$109M⊙ and a massive building block
of mass M${}_{\star}=$ 1.43$\times$109M⊙). As can be seen in Fig 11, there are
subtle differences in the spatial distribution (at $z=0$) for all populations
examined between m12i and m12m. More specifically, m12i presents a clear thin
young disc population in addition to a spherical old disc, main branch, and
building block (of mass M${}_{\star}=$ 2.67$\times$107M⊙). Conversely, m12m
presents a much flatter distribution, where the main branch and massive
building block systems present an oval shape, in a similar fashion to the old
disc in this Milky Way-mass galaxy.

Moreover, we find that in m12i, the main branch and building block overlap in
the Toomre diagram and are largely non-rotational ($v_{\phi}$$\sim$0 km/s),
and are different to the old/young disc, which show a more extended
distribution rotating at higher $v_{\phi}$. Conversely, m12m reveals that its
main branch and massive building block also overlap in kinematic space, but
rotate at a higher tangential velocity of $v_{\phi}$$\sim$100 km/s, lagging
behind the young disc by $\sim$150 km/s. In this halo, the main branch and
massive building block also overlap with the old disc. The difference in the
$v_{\phi}$ magnitude between the proto-Milky Way system in m12i and m12m is
interesting. However, it is likely due to m12i being an outlier in terms of
the rotational velocities of its metal-poor stars (Santistevan et al., 2022).
All the other Latte Milky Way-mass galaxies show a preference for prograde
disc orbits for older and/or lower metallicity stars, similar to m12m
(although this system also shows higher than average rotation, see the right
panel of Fig 7). Santistevan et al. (2022) argue that metal-poor/old stars on
prograde disc orbits are a consequence of major building blocks depositing
stars/gas on prograde orbits, which typically set the orientation of the
resulting Milky Way-mass galaxy disc at $z=0$.

Figure 12: Age metallicity relation for the same samples as shown in Fig 11.
For this particular comparison we see that m12i, a Milky Way-mass halo with a
proto-galaxy population formed primarily by one massive and dominant halo
(i.e, main branch), is older than the proto-galaxy population of m12m, a Milky
Way-mass halo formed from two dominant halos (i.e., a main branch and a
massive building block). We argue that this $\sim 2$ Gyr difference should be
detectable in relative ages of metal-poor stars. Moreover, the peak in the
metallicity distribution for the proto-Milky Way population in m12i is more
metal-poor than the one in m12m, on the order of $\sim$1 dex. These subtle
differences are also seen for $\sim 60\%$ of the parent Milky Way-mass galaxy
sample (see Section 5.5 for details). These results hint that age estimates
and [Fe/H] information for large samples of metal-poor stars in the inner
regions of the Milky Way could possibly help answer if the Galaxy formed from
one dominant halo or from two.

In terms of their [Mg/Fe]-[Fe/H] compositions, we find again that for both
m12i and m12m, the main branch and building block populations overlap. This is
to be expected for m12m (see Fig 8 and Fig 9), but is a surprise for m12i.
Interestingly, the [Mg/Fe]-[Fe/H] compositions of these systems are clearly
different from that of the old/young discs in m12i, which present a much
higher [Fe/H] distribution. However, we note that there is an overlap between
the metal-rich sequence of the main branch/building block and the metal-poor
tail of the old disc. Conversely, for m12m this overlap between the main
branch and massive building block is more pronounced. This highlights that the
dominant proto-Milky Way populations are likely to overlap with the metal-
poor/old disc (Horta et al., 2021a; Conroy et al., 2021; Mardini et al., 2022;
Horta et al., 2023b).

In a similar vein, another interesting diagram to investigate relates to the
time difference between proto-Milky Way systems that form from one main system
versus two main haloes. To investigate this point further, in Fig 12 we show
the main branch, (massive) building block, and old/young disc in the age-
metallicity plane. Here the data is displayed in the same way as in Fig 11. We
find that for this particular comparison, a Milky Way-mass halo with a proto-
galaxy population formed primarily by one massive and dominant main branch
halo, m12i (but also Thelma, Louise, Romeo, and Juliet, see Fig 3) is older
than the proto-galaxy population of a Milky Way-mass halo formed from one
dominant main branch halo and a massive building block halo, m12m (but also
m12c and Romulus). This difference is on the order of $\sim$2 Gyr.
Furthermore, the peak in the metallicity distribution for the proto-Milky Way
population in m12i is more metal-poor than the one in m12m, on the order of
$\sim$1 dex. Although this comparison has only been shown for two haloes, we
find that it is qualitatively satisfied for 8 of the 13 Milky Way-mass
galaxies. We argue that this result is potentially very important, and
suggests that age differences in metal-poor stars in the inner regions of
Milky Way-mass haloes could possibly help decipher if proto-Milky Way systems
formed from one or two main haloes. However, we do note that for m12b, m12f,
m12r, m12w, and Remus, this is not the case.

In summary, we argue that disentangling if a proto-Milky Way system formed
from one or two dominant systems is going to be challenging, as their
obervable properties are going to be similar. However, we suggest that with
large and complete samples of metal-poor stars, that include elemental
abundances synthesised in nucleosynthetic sites not explored in this work, as
well as accurate age estimates, may hold the clues to answering this question.

###  5.6 Hunting for the proto-galaxy in the Milky Way

Given all our findings, we now provide a list of ideas/pointers we suggest one
to follow if aiming to identify the majority of the populations belonging to
the proto-Milky Way:

i) Mass and age: we suggest that the best possible way to find the proto-
Galaxy in the Milky Way would be to identify stellar populations older than
$\tau$$\gtrsim$10 Gyr, that amount to a mass smaller than
M⋆$\lesssim$109–1010M⊙.

ii) Spatial distribution: this stellar population would primarily be
concentrated towards the inner regions of the Galaxy (95$\%$ of the mass
within $r\sim 30-40$ kpc), with 50$\%$ of the mass contained within $\sim
5-10$ kpc. It would also likely have an oblate profile
($\varepsilon$$\sim$0.7).

iii) Kinematics: the proto-Milky Way would not be a strongly rotating
population (w.r.t the present day Galactic disc,
0.25$\lesssim$$\langle$$\kappa$/$\kappa_{\mathrm{disc}}$$\rangle$
$\lesssim$0.8). However, it is still likely to host prograde orbits and
manifest moderate tangential velocity values
(0$\lesssim$$v_{\phi}$$\lesssim$60 km/s, see Fig 7), matching that of an old
and metal-poor disc (Fig 11).

iv) Chemical compositions: the proto-Milky Way could have a wide range of
[Fe/H] and [$\alpha$/Fe] values. We find that the majority of the stars
comprising the proto-Milky Way (either from one main branch system or a main
branch and a massive building block) are likely to overlap with the metal-weak
old disc in [$\alpha$/Fe] and [Fe/H]. However, additional elemental abundance
ratios, which have not been explored in this work, may provide additional
separation between proto-Galaxy fragments and other co-spatial populations.

Hunts for the proto-Galaxy and/or the building blocks that formed it have been
performed recently in observational studies (Horta et al., 2021a; Belokurov &
Kravtsov, 2022; Rix et al., 2022) using chemical-kinematic information. These
results are shedding light on the earliest stages of formation of the Milky
Way. However, open questions still remain, such as: how many systems comprise
the proto-Milky Way? when did the proto-Milky Way emerge? what caused the
proto-Galaxy to emerge and morph into the metal-poor disc? We believe that
this study has provided some intuition on how to answer some of these
questions, as we have been able to assess: a) what proto-Milky Way populations
are made of; b) how massive and/or old proto-Milky Way systems are; c) how and
if different proto-Milky Way fragments can be distinguishable in chemical-
kinematic samples.

The advent of the SDSS-V project Milky Way Mapper (Kollmeier et al., 2017)
will deliver precise elemental abundance ratios for over $\sim$50,000 stars
with [Fe/H]$<-$1, that when complemented with the Gaia mission, and other
upcoming massive spectroscopic surveys (e.g., WEAVE: Dalton et al., 2012;
4MOST: de Jong et al., 2019), will provide an unprecedented amount of
chemical-kinematic information for metal-poor stars in the innermost Galaxy.
Along those lines, recent work with the Gaia XP spectra are delivering a
colossal amount of metallicities for stars in the inner regions of the Galaxy
(Andrae et al., 2023), which are also helping resolve populations in this
region (Rix et al., 2022). In a similar vein, detailed spectroscopic follow up
(e.g., PRISTINE: Starkenburg et al. 2017) of metal-poor stars in the innermost
regions of the Galaxy can also help disentangle the earliest stages of
formation of the Milky Way (Arentsen et al., 2020b, a; Lucey et al., 2020;
Sestito et al., 2023).

Given our results, we suggest that with these data it should be possible to
pick out the massive contributors to proto-Milky Way systems from the smaller
building blocks. We are only scratching the surface, but we are on the right
path to finding the Milky Way’s heart and understanding its early assembly
history.

##  6 Limitations and further considerations

Before listing our conclusions, it is imperative we discuss the limitations of
our work. Below we provide a summary of the limitations and considerations for
the interpretations of our findings:

  * •

The new aspect of this work has been the ability to track the systems that
comprise proto-Milky Way populations. However, in doing so, we have decided to
track only those luminous halos that are resolvable given the resolution
limits of FIRE-2. Specifically, we have tracked the populations from halos
with 150 star particles or more, leading to a minimum stellar mass of
$\sim$1$\times$106 M⊙ for the Latte suite, and $\sim$5$\times$105 M⊙ in the
ELVIS suite101010The small difference in the resolution between the two suites
is small enough to not have any impact on our results, as we are tracking
haloes with 150 star particles or more.. Albeit our hands being tied with the
ability to resolve haloes, we have been cautious in our choices to be able to
track as many resolvable systems as possible in these simulations. However, it
is likely that in reality there are more luminous halos of lower masses that
contribute to the build-up of proto-Milky Way populations. We argue that the
contribution in mass of these lower mass building blocks is small, but cannot
rule out the possibility of it being non-negligible (see Santistevan et al.,
2020 and Ghandi et al ., in prep. for an accounting of lower mass building
blocks). Furthermore, the kinematic properties of the lower mass haloes may
not be fully resolved with 150 star particles. Although this only affects a
small fraction of our sample, it is another limitation to keep in mind.
Studies focusing on tackling the assembly history of proto-Milky Way
populations, and the properties of their constituent main branch and building
block systems, may need to take this issue into account when comparing to our
findings.

  * •

A key property we have defined in this work has been $t_{\mathrm{MR_{3:1}}}$,
the time at which a proto-Milky Way emerges. Although well reasoned (see
Section 3 for details), this choice is arbitrary, and has ramifications on all
the present day properties of the proto-Milky Way populations, as well as the
number of building block systems each proto-Milky Way inherits. It could have
also been equally as valid to define this time in another well motivated way
(for example, the time in which the host halo peaks in star formation). This
is beyond the scope of this work, but could be an interesting way of expanding
on these results.

  * •

Here we have only studied the properties of the field stellar components of
proto-Milky Way populations. However, galaxies across the Cosmos contain
globular clusters (GCs). In the Milky Way, it has been shown that the
disruption of GCs contribute a significant amount to the total stellar halo
mass budget (Martell et al., 2017; Schiavon et al., 2017; Koch et al., 2019;
Horta et al., 2021b; Belokurov & Kravtsov, 2023). This would imply that in
addition to the field components comprising proto-Milky Way systems, one must
also take into account the contribution from disrupted/evaporated GCs. For
this work, this limitation is due to FIRE-2 not including prescriptions for
the formation and evolution of GCs. However, we argue that it is an important
point that needs to be considered for observational and potential future
simulation work.

  * •

When comparing the observable properties of main branch and building block
events, we have primarily only examined qualitatively the mean/median values.
It would be interesting to investigate how the average and full distribution
of the spatial, kinematic, and chemical properties of these systems compare
quantitatively (following methods and tools that already exist; e.g.,
Cunningham et al., 2022; Horta et al., 2023a).

##  7 Conclusions

At the earliest stages of formation, galaxies experience rapid and chaotic
growth, either by coalescence of low mass galaxies/clumps and/or filamentary
supply of gas. In the FIRE-2 simulations, bursty stellar feedback that
repeatedly blows apart the ISM at early times also appears to be critical to
set the properties of early stellar populations (Yu et al., 2021; Gurvich et
al., 2023; Hopkins et al., 2023, e.g.,). The remains of the stars born during
this phase constitute the proto-galaxy, and should retain the clues to
understanding galaxy formation at these earliest stages. In this work, we have
searched for the fragments (namely, the main branch and building blocks) that
constitute the proto-galaxy in thirteen Milky Way-mass haloes from the FIRE-2
simulations. We then examined their observable properties at present day, with
the aim of answering the following questions:

  * •

What constitutes a proto-Milky Way? We find that proto-Milky Way populations
are made of either one ($\sim 60\%$) or two ($\sim 40\%$) dominant systems of
similar mass to the LMC (i.e., M⋆$\sim$1$\times$109M⊙) and $\sim$3-5 other
lower-mass building blocks with an average stellar mass of
M⋆$\sim$4$\times$107M⊙ (see Fig 2, Fig 3, and Fig 10). The case of two clear
dominant systems comprising the proto-Milky Way is especially clear in m12f
and m12m. However, the number of building blocks we find in this study that
comprise a proto-Milky Way is grounded by our choice to only track systems
with 150 star particles or more.

  * •

When does a proto-Milky Way form? Given our assumptions (Section 3) we find
that on average, proto-Milky Way populations are old ($\tau$$\gtrsim$10 Gyr or
$z$$\gtrsim$2, see Fig 3).

  * •

Does environment play a role? The noticeable differences found between
galaxies in different environments are the times in which proto-Milky Way
systems assemble, and the sizes of the average building block. Table 2 shows
that on average, systems in pairs assemble earlier and from smaller mass
systems than proto-Milky Way populations in isolation, in line with results
from Santistevan et al. (2020). We also find that pairs tend to contain 50$\%$
of their stellar populations closer to the host’s centre when compared to
isolated systems (Fig 5), and that pairs tend to present slightly more
enriched median [$\alpha$/Fe] and/or [Fe/H] composition values.

  * •

Where are the debris of proto-Milky Way systems spatially contained? The
dominant components of the proto-Milky Way (namely, the main branch population
and the most massive building blocks) contain 50$\%$ of their stellar mass
within $r\lesssim 5-10$ kpc, and 95$\%$ within $r\lesssim 40-60$ kpc (see Fig
4 and Fig 5).

  * •

What shape do debris from the proto-Milky Way have? Although different main
branch/building blocks adopt a range of morphologies, the dominant components
of proto-Milky Way systems adopt an oblate shape ($\varepsilon$$\sim$0.7, see
Fig 6).

  * •

What kinematics do the debris of proto-Milky Way systems have? All the
fragments of proto-Milky Way systems show low level of systematic net rotation
with respect to the present day disc of the host galaxy (i.e.,
$\langle\kappa$/$\kappa_{\mathrm{disc}}\rangle\lesssim 0.8$), but are also not
purely isotropic (Fig 7). The majority of these stellar populations rotate on
prograde orbits, and can reach average azimuthal velocities of up to
$v_{\phi}\sim 100-150$ km/s.

  * •

What are the chemical compositions of proto-Milky Ways? The main branch and
building block components of proto-Milky Way systems can present a wide range
of [Fe/H] and [$\alpha$/Fe] values, and can adopt high [$\alpha$/Fe] values
(Fig 8). This would make it difficult to distinguish the dominant proto-Milky
Way populations from populations formed later in the disc using the average
value of these abundances alone (see Fig 11). However, disentangling the
massive fragments from the lower mass ones is possible with [$\alpha$/Fe] and
[Fe/H] compositions.

  * •

Can we distinguish the fragments that build up proto-Milky Way systems? On the
whole, it is possible to separate the dominant components of proto-Milky Way
systems (i.e., the main branch) from the non-dominant (low mass building
blocks) components using chemical-kinematic samples. However, distinguishing
stars formed in the main branch progenitor from the most massive building
blocks will likely be difficult owing to the big overlap in all chemical and
kinematic planes shown in this work. However, studies focusing on other
chemical abundance measurements and/or studies with larger samples may be able
to disentangle the dominant proto-Milky Way fragments based on chemical
abundance measurements.

  * •

How do we hunt for the proto-Milky Way? We suggest that metal-poor stars
confined to the inner Galaxy are likely to be the easiest to find, as they are
the debris from the more dominant fragments that formed the proto-Milky Way.
These stars should have prograde orbits that are not strongly rotating with
the Milky Way’s disc, but are also not purely isotropic. They should present,
on average, more enriched chemical abundance ratios when compared to lower
mass building blocks and later accreted satellite galaxies, and are likely to
overlap in chemical space with the metal-poor tail of the old disc. Thus, we
suggest that additional chemical abundance information, likely probing
different nucleosynthetic channels may help disentangle different fragments of
the proto-Milky Way.

  * •

Is the proto-Milky Way formed from one or two dominant haloes? It may be
possible to answer this question by examining relative age differences between
metal-poor stars in the central regions of the Galaxy. From our comparison of
m12i and m12m in Fig 12, proto-Milky Way systems formed from one dominant halo
may tend to assemble earlier, and are thus older and more metal-poor, when
compared to proto-Milky Way systems formed from two dominant haloes. However,
we do stress that this result is from one comparison alone, and is only
applicable to $\sim 60\%$ of the full sample (8 of the 13 Milky Way-mass
galaxies studied).

## Acknowledgements

DH thanks Ricardo P. Schiavon, Vasily Belokurov, Hans-Walter Rix, and
Stephanie Monty for helpful discussions. He also thanks Sue, Alex, and Debra
for their constant support. ECC acknowledges support for this work provided by
NASA through the NASA Hubble Fellowship Program grant HST-HF2-51502 awarded by
the Space Telescope Science Institute, which is operated by the Association of
Universities for Research in Astronomy, Inc., for NASA, under contract
NAS5-26555. RES gratefully acknowledges support from NSF grant AST-2007232 and
NASA grant 19-ATP19-0068. CAFG was supported by NSF through grants AST-2108230
and CAREER award AST-1652522; by NASA through grants 17-ATP17-0067 and
21-ATP21-0036; by STScI through grant HST-GO-16730.016-A; and by CXO through
grant TM2-23005X. PJG received support from the Texas Advanced Computing
Center (TACC) via a Frontera Computational Science Fellowship.

## Data Availability

FIRE-2 simulations are publicly available (Wetzel et al., 2023) at
http://flathub.flatironinstitute.org/fire. Additional FIRE simulation data is
available at https://fire.northwestern.edu/data. A public version of the Gizmo
code is available at http://www.tapir.caltech.edu/~phopkins/Site/GIZMO.html.

## References

  * Abadi et al. (2006) Abadi M. G., Navarro J. F., Steinmetz M., 2006, MNRAS, 365, 747
  * Amorisco (2017) Amorisco N. C., 2017, MNRAS, 464, 2882
  * Andrae et al. (2023) Andrae R., Rix H.-W., Chandra V., 2023, arXiv e-prints, p. arXiv:2302.02611
  * Arentsen et al. (2020a) Arentsen A., et al., 2020a, MNRAS, 491, L11
  * Arentsen et al. (2020b) Arentsen A., et al., 2020b, MNRAS, 496, 4964
  * Behroozi et al. (2013) Behroozi P. S., Wechsler R. H., Wu H.-Y., 2013, ApJ, 762, 109
  * Bell et al. (2008) Bell E. F., et al., 2008, ApJ, 680, 295
  * Bellardini et al. (2021) Bellardini M. A., Wetzel A., Loebman S. R., Faucher-Giguère C.-A., Ma X., Feldmann R., 2021, MNRAS, 505, 4586
  * Belokurov & Kravtsov (2022) Belokurov V., Kravtsov A., 2022, MNRAS, 514, 689
  * Belokurov & Kravtsov (2023) Belokurov V., Kravtsov A., 2023, arXiv e-prints, p. arXiv:2306.00060
  * Belokurov et al. (2007) Belokurov V., et al., 2007, ApJ, 658, 337
  * Belokurov et al. (2018) Belokurov V., Erkal D., Evans N. W., Koposov S. E., Deason A. J., 2018, MNRAS, 478, 611
  * Belokurov et al. (2020) Belokurov V., Sanders J. L., Fattahi A., Smith M. C., Deason A. J., Evans N. W., Grand R. J. J., 2020, MNRAS, 494, 3880
  * Bonaca et al. (2017) Bonaca A., Conroy C., Wetzel A., Hopkins P. F., Kereš D., 2017, ApJ, 845, 101
  * Brown et al. (2006) Brown T. M., Smith E., Ferguson H. C., Rich R. M., Guhathakurta P., Renzini A., Sweigart A. V., Kimble R. A., 2006, ApJ, 652, 323
  * Bullock & Johnston (2005) Bullock J. S., Johnston K. V., 2005, ApJ, 635, 931
  * Carollo et al. (2007) Carollo D., et al., 2007, Nature, 450, 1020
  * Chiba & Beers (2000) Chiba M., Beers T. C., 2000, AJ, 119, 2843
  * Conroy et al. (2021) Conroy C., Naidu R. P., Garavito-Camargo N., Besla G., Zaritsky D., Bonaca A., Johnson B. D., 2021, Nature, 592, 534
  * Cooper et al. (2010) Cooper A. P., et al., 2010, MNRAS, 406, 744
  * Cunningham et al. (2022) Cunningham E. C., et al., 2022, ApJ, 934, 172
  * Dalton et al. (2012) Dalton G., et al., 2012, in McLean I. S., Ramsay S. K., Takami H., eds, Society of Photo-Optical Instrumentation Engineers (SPIE) Conference Series Vol. 8446, Ground-based and Airborne Instrumentation for Astronomy IV. p. 84460P, doi:10.1117/12.925950
  * Deason et al. (2011) Deason A. J., Belokurov V., Evans N. W., 2011, MNRAS, 416, 2903
  * Deason et al. (2023) Deason A. J., Koposov S. E., Fattahi A., Grand R. J. J., 2023, MNRAS, 520, 6091
  * Di Matteo et al. (2019) Di Matteo P., Haywood M., Lehnert M. D., Katz D., Khoperskov S., Snaith O. N., Gómez A., Robichon N., 2019, A&A, 632, A4
  * Donlon et al. (2022) Donlon Thomas I., Newberg H. J., Kim B., Lépine S., 2022, ApJ, 932, L16
  * El-Badry et al. (2018) El-Badry K., et al., 2018, MNRAS, 480, 652
  * Escala et al. (2018) Escala I., et al., 2018, MNRAS, 474, 2194
  * Escala et al. (2019) Escala I., Kirby E. N., Gilbert K. M., Cunningham E. C., Wojno J., 2019, ApJ, 878, 42
  * Escala et al. (2020) Escala I., Gilbert K. M., Kirby E. N., Wojno J., Cunningham E. C., Guhathakurta P., 2020, ApJ, 889, 177
  * Faucher-Giguère et al. (2009) Faucher-Giguère C.-A., Lidz A., Zaldarriaga M., Hernquist L., 2009, ApJ, 703, 1416
  * Ferguson et al. (2002) Ferguson A. M. N., Irwin M. J., Ibata R. A., Lewis G. F., Tanvir N. R., 2002, AJ, 124, 1452
  * Font et al. (2008) Font A. S., Johnston K. V., Ferguson A. M. N., Bullock J. S., Robertson B. E., Tumlinson J., Guhathakurta P., 2008, ApJ, 673, 215
  * Font et al. (2011) Font A. S., McCarthy I. G., Crain R. A., Theuns T., Schaye J., Wiersma R. P. C., Dalla Vecchia C., 2011, MNRAS, 416, 2802
  * Forbes (2020) Forbes D. A., 2020, MNRAS, 493, 847
  * Fragkoudi et al. (2020) Fragkoudi F., et al., 2020, MNRAS, 494, 5936
  * Gaia Collaboration et al. (2022) Gaia Collaboration et al., 2022, arXiv e-prints, p. arXiv:2208.00211
  * Gandhi et al. (2022) Gandhi P. J., Wetzel A., Hopkins P. F., Shappee B. J., Wheeler C., Faucher-Giguère C.-A., 2022, arXiv e-prints, p. arXiv:2202.10477
  * Garrison-Kimmel et al. (2017) Garrison-Kimmel S., et al., 2017, MNRAS, 471, 1709
  * Garrison-Kimmel et al. (2019a) Garrison-Kimmel S., et al., 2019a, MNRAS, 487, 1380
  * Garrison-Kimmel et al. (2019b) Garrison-Kimmel S., et al., 2019b, MNRAS, 489, 4574
  * Grillmair & Dionatos (2006) Grillmair C. J., Dionatos O., 2006, ApJ, 643, L17
  * Gurvich et al. (2023) Gurvich A. B., et al., 2023, MNRAS, 519, 2598
  * Hahn & Abel (2011) Hahn O., Abel T., 2011, MNRAS, 415, 2101
  * Haywood et al. (2018) Haywood M., Di Matteo P., Lehnert M. D., Snaith O., Khoperskov S., Gómez A., 2018, ApJ, 863, 113
  * Helmi (2020) Helmi A., 2020, arXiv e-prints, p. arXiv:2002.04340
  * Helmi et al. (2018) Helmi A., Babusiaux C., Koppelman H. H., Massari D., Veljanoski J., Brown A. G. A., 2018, Nature, 563, 85
  * Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
  * Hopkins (2016) Hopkins P. F., 2016, MNRAS, 455, 89
  * Hopkins et al. (2018a) Hopkins P. F., et al., 2018a, MNRAS, 477, 1578
  * Hopkins et al. (2018b) Hopkins P. F., et al., 2018b, MNRAS, 480, 800
  * Hopkins et al. (2023) Hopkins P. F., et al., 2023, MNRAS,
  * Horta et al. (2021a) Horta D., et al., 2021a, MNRAS, 500, 1385
  * Horta et al. (2021b) Horta D., et al., 2021b, MNRAS, 500, 5462
  * Horta et al. (2023a) Horta D., et al., 2023a, MNRAS, 520, 5671
  * Horta et al. (2023b) Horta D., et al., 2023b, ApJ, 943, 158
  * Ibata et al. (1994) Ibata R. A., Gilmore G., Irwin M. J., 1994, Nature, 370, 194
  * Johnston et al. (2008) Johnston K. V., Bullock J. S., Sharma S., Font A., Robertson B. E., Leitner S. N., 2008, ApJ, 689, 936
  * Khoperskov et al. (2022a) Khoperskov S., et al., 2022a, arXiv e-prints, p. arXiv:2206.04521
  * Khoperskov et al. (2022b) Khoperskov S., et al., 2022b, arXiv e-prints, p. arXiv:2206.04522
  * Kirby et al. (2013) Kirby E. N., Cohen J. G., Guhathakurta P., Cheng L., Bullock J. S., Gallazzi A., 2013, ApJ, 779, 102
  * Kirby et al. (2020) Kirby E. N., Gilbert K. M., Escala I., Wojno J., Guhathakurta P., Majewski S. R., Beaton R. L., 2020, AJ, 159, 46
  * Koch et al. (2019) Koch A., Grebel E. K., Martell S. L., 2019, A&A, 625, A75
  * Kollmeier et al. (2017) Kollmeier J. A., et al., 2017, SDSS-V: Pioneering Panoptic Spectroscopy (arXiv:1711.03234)
  * Kroupa (2001) Kroupa P., 2001, MNRAS, 322, 231
  * Kruijssen et al. (2020) Kruijssen J. M. D., et al., 2020, arXiv e-prints, p. arXiv:2003.01119
  * Krumholz & Gnedin (2011) Krumholz M. R., Gnedin N. Y., 2011, ApJ, 729, 36
  * Leitherer et al. (1999) Leitherer C., et al., 1999, ApJS, 123, 3
  * Lucey et al. (2020) Lucey M., et al., 2020, The COMBS Survey II: Distinguishing the Metal-Poor Bulge from the Halo Interlopers (arXiv:2009.03886)
  * Ma et al. (2017) Ma X., Hopkins P. F., Wetzel A. R., Kirby E. N., Anglés-Alcázar D., Faucher-Giguère C.-A., Kereš D., Quataert E., 2017, MNRAS, 467, 2430
  * Mackereth et al. (2019) Mackereth J. T., et al., 2019, MNRAS, 482, 3426
  * Mardini et al. (2022) Mardini M. K., Frebel A., Chiti A., Meiron Y., Brauer K. V., Ou X., 2022, ApJ, 936, 78
  * Martell et al. (2017) Martell S. L., et al., 2017, MNRAS, 465, 3203
  * McCarthy et al. (2012) McCarthy I. G., Font A. S., Crain R. A., Deason A. J., Schaye J., Theuns T., 2012, MNRAS, 420, 2245
  * McCluskey et al. (2023) McCluskey F., Wetzel A., Loebman S. R., Moreno J., Faucher-Giguere C.-A., 2023, arXiv e-prints, p. arXiv:2303.14210
  * Necib et al. (2019) Necib L., Lisanti M., Garrison-Kimmel S., Wetzel A., Sanderson R., Hopkins P. F., Faucher-Giguère C.-A., Kereš D., 2019, ApJ, 883, 27
  * Newberg et al. (2009) Newberg H. J., Yanny B., Willett B. A., 2009, ApJ, 700, L61
  * Orkney et al. (2023) Orkney M. D. A., et al., 2023, arXiv e-prints, p. arXiv:2303.02147
  * Panithanpaisal et al. (2021) Panithanpaisal N., Sanderson R. E., Wetzel A., Cunningham E. C., Bailin J., Faucher-Giguère C.-A., 2021, arXiv e-prints, p. arXiv:2104.09660
  * Planck Collaboration et al. (2020) Planck Collaboration et al., 2020, A&A, 641, A6
  * Rix et al. (2022) Rix H.-W., et al., 2022, ApJ, 941, 45
  * Samuel et al. (2020) Samuel J., et al., 2020, MNRAS, 491, 1471
  * Sanderson et al. (2018) Sanderson R. E., et al., 2018, ApJ, 869, 12
  * Sanderson et al. (2020) Sanderson R. E., et al., 2020, ApJS, 246, 6
  * Santistevan et al. (2020) Santistevan I. B., Wetzel A., El-Badry K., Bland-Hawthorn J., Boylan-Kolchin M., Bailin J., Faucher-Giguère C.-A., Benincasa S., 2020, MNRAS, 497, 747
  * Santistevan et al. (2022) Santistevan I. B., Wetzel A., Tollerud E., Sanderson R. E., Samuel J., 2022, arXiv e-prints, p. arXiv:2208.05977
  * Schiavon et al. (2017) Schiavon R. P., et al., 2017, MNRAS, 465, 501
  * Semenov et al. (2023) Semenov V. A., Conroy C., Chandra V., Hernquist L., Nelson D., 2023, arXiv e-prints, p. arXiv:2306.13125
  * Sestito et al. (2023) Sestito F., et al., 2023, MNRAS, 518, 4557
  * Shipp et al. (2022) Shipp N., et al., 2022, arXiv e-prints, p. arXiv:2208.02255
  * Starkenburg et al. (2017) Starkenburg E., et al., 2017, MNRAS, 471, 2587
  * Su et al. (2017) Su K.-Y., Hopkins P. F., Hayward C. C., Faucher-Giguère C.-A., Kereš D., Ma X., Robles V. H., 2017, MNRAS, 471, 144
  * Wetzel & Garrison-Kimmel (2020a) Wetzel A., Garrison-Kimmel S., 2020a, HaloAnalysis: Read and analyze halo catalogs and merger trees (ascl:2002.014)
  * Wetzel & Garrison-Kimmel (2020b) Wetzel A., Garrison-Kimmel S., 2020b, GizmoAnalysis: Read and analyze Gizmo simulations (ascl:2002.015)
  * Wetzel et al. (2016) Wetzel A. R., Hopkins P. F., Kim J.-h., Faucher-Giguère C.-A., Kereš D., Quataert E., 2016, ApJ, 827, L23
  * Wetzel et al. (2023) Wetzel A., et al., 2023, ApJS, 265, 44
  * Wheeler et al. (2019) Wheeler C., et al., 2019, MNRAS, 490, 4447
  * White & Frenk (1991) White S. D. M., Frenk C. S., 1991, ApJ, 379, 52
  * White & Rees (1978) White S. D. M., Rees M. J., 1978, MNRAS, 183, 341
  * Yu et al. (2021) Yu S., et al., 2021, MNRAS, 505, 889
  * Yu et al. (2023) Yu S., et al., 2023, MNRAS, 523, 6220
  * de Jong et al. (2019) de Jong R. S., et al., 2019, The Messenger, 175, 3