# The history of the dark and luminous side of Milky Way-like progenitors

Luca Graziani, M. de Bennassuti, R. Schneider, D. Kawata, S. Salvadori

► To cite this version:

Luca Graziani, M. de Bennassuti, R. Schneider, D. Kawata, S. Salvadori. The history of the dark and luminous side of Milky Way-like progenitors. Monthly Notices of the Royal Astronomical Society, 2017, 469 (1), pp.1101-1116. 10.1093/mnras/stx900. hal-01582453

HAL Id: hal-01582453

https://hal.science/hal-01582453v1

Submitted on 29 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
# The history of the dark and luminous side of Milky Way-like progenitors

L. Graziani,<fnref n="1" />★ M. de Bennassuti,<fnref n="1" />,<fnref n="2" /> R. Schneider,<fnref n="1" />,<fnref n="2" /> D. Kawata<fnref n="3" /> and S. Salvadori<fnref n="4" />

Accepted 2017 April 10. Received 2017 April 9; in original form 2016 December 4

## ABSTRACT

Here we investigate the evolution of a Milky Way (MW)-like galaxy with the aim of predicting the properties of its progenitors all the way from $z \sim 20$ to $z = 0$. We apply GAMESH to a high-resolution N-body simulation following the formation of a MW-type halo and we investigate its properties at $z \sim 0$ and its progenitors in $0 < z < 4$. Our model predicts the observed galaxy main sequence, the mass-metallicity and the Fundamental Plane of metallicity relations in $0 < z < 4$. It also reproduces the stellar mass evolution of candidate MW progenitors in $0 \lesssim z \lesssim 2.5$, although the star formation rate and gas fraction of the simulated galaxies follow a shallower redshift dependence. We find that while the MW star formation and chemical enrichment are dominated by the contribution of galaxies hosted in Lyman $\alpha$ cooling haloes, at $z > 6$ the contribution of star-forming minihaloes is comparable to the star formation rate along the MW merger tree. These systems might then provide an important contribution in the early phases of reionization. A large number of minihaloes with old stellar populations, possibly Population III stars, are dragged into the MW or survive in the Local Group. At low redshift dynamical effects, such as halo mergers, tidal stripping and halo disruption redistribute the baryonic properties among halo families. These results are critically discussed in light of future improvements including a more sophisticated treatment of radiative feedback and inhomogeneous metal enrichment.

**Key words:** stars: formation – stars: Population II – galaxies: evolution – galaxies: formation – galaxies: stellar content – dark ages, reionization, first stars – cosmology: theory.

## 1 INTRODUCTION

Modern cosmological models (White & Frenk 1991; Cole et al. 1994, 2000; Henriques et al. 2015) interpret the properties of galaxies observed in the present Universe as the result of the intricate interplay in feedback mechanisms acting during halo mass assembly and shaping the galactic baryons through cosmic times (Mo, van den Bosch & White 2010). Despite the increasing number of objects provided by large-scale surveys (Huchra et al. 1983; Abazajian et al. 2003; Alam et al. 2015; Blanton et al. 2017) and by recent high-redshift observations (Grogin et al. 2011; Bouwens et al. 2016), the incomparable level of details available at small scale in our Galaxy still places the Milky Way (MW) at the centre of any model of galaxy formation and evolution.

The possibility of resolving stars, both in the MW and in the closest galaxies of the observed Local Group (oLG), provides a unique observational data set and allows to build Galactic archaeology

models on solid observational grounds (Scannapieco et al. 2006; Tumlinson 2006, 2010; Brook et al. 2007; Komiya et al. 2007; Salvadori et al. 2010a; Frebel & Bromm 2012; de Bennassuti et al. 2014; Hartwig et al. 2015; Salvadori, Skúladóttir & Tolstoy 2015; de Bennassuti et al. 2017).

The great number of kinematic and chemical tracers available in the MW (see Bland-Hawthorn & Gerhard 2016 for a recent review) has been complemented by detailed observations of stellar populations in nearby dwarf galaxies, enabling us to infer their star formation histories and to interpret the nature of the smallest objects (Tolstoy, Hill & Tosi 2009; McConnellie 2012; Madden et al. 2013; Monelli et al. 2016). Ultrafaint dwarf galaxies (Simon & Geha 2007; Kirby et al. 2008; Brown et al. 2014), for example, are believed to be fossil remnants of the pre-reionization era and the record of radiative feedback at play in the early Universe (Bovill & Ricotti 2009; Muñoz et al. 2009; Salvadori & Ferrara 2009; Vincenzo et al. 2014; Bland-Hawthorn, Sutherland & Webster 2015; Salvadori et al. 2015). Given this wealth of data, theoretical models have attempted to understand how the star formation history of the MW and its dwarf companions (Kennicutt & Evans 2012) is

<footnote n="1">¹INAF-Osservatorio Astronomico di Roma, Via di Frascati 33, I-00078 Monte Porzio Catone, Italy</footnote>

<footnote n="2">²Dipartimento di Fisica, Sapienza, Università di Roma, Piazzale Aldo Moro 5, I-00185 Roma, Italy</footnote>

<footnote n="3">³Mullard Space Science Laboratory, University College London, Holmbury St. Mary, Dorking, Surrey RH5 6NT, UK</footnote>

<footnote n="4">⁴GEPI, Observatoire de Paris, PSL Research University, CNRS, Univ Paris Diderot, Sorbonne Paris Cité, Place Jules Janssen, F-92195 Meudon, France</footnote>

<footnote>*E-mail: luca.graziani@oa-roma.inaf.it</footnote>
affected by large-scale processes like cosmic reionization and metal enrichment (Ocvirk et al. 2014; Graziani et al. 2015; Wetzel, Deason & Garrison-Kimmel 2015).

Self-consistent models across scales and cosmic times tuned on low-redshift observations have the enormous potential of placing stringent constraints on the nature of unobserved objects at high redshift, as for example the first generation of stars and galaxies (Bromm & Yoshida 2011), also helping to investigate the efficiency of early radiative and chemical processes and their mutual impact (Ciardi & Ferrara 2005).

Large-scale dark matter (DM) simulations, as the Millennium suite<fnref n="1" /> (Springel et al. 2005; Springel, Frenk & White 2006; Boylan-Kolchin et al. 2009; Angulo et al. 2012) or the MultiDark and Bolshoi runs<fnref n="2" /> (Klypin, Trujillo-Gomez & Primack 2011), have traditionally provided the theoretical framework on top of which semi-analytic models (SAMs) interpret global observed quantities (see for instance Knebe et al. 2015). Recently, two hydrodynamical projects (Vogelsberger et al. 2014; Schaye et al. 2015) including detailed feedback on subgrid scales have been able to reach an unprecedented realism in reproducing morphological and structural galaxy properties.

In the past years many DM simulations of MW-like haloes and of the oLG have been successfully performed, often as spin-off of large-scale simulations, such as the *AQUARIUS* run (Springel et al. 2008), the *VIA LACTEA* project (Diemand et al. 2008), the *CLUES* project (Nuza et al. 2014; Yepes, Gottlöber & Hoffman 2014) and the *ELVIS* simulation suite (Garrison-Kimmel et al. 2014). They have been extensively used to study the structural properties of DM haloes, the statistics of their satellites, as well as to correctly constrain the initial conditions leading to the dynamical configuration of our local Universe, having an MW–M31 galaxy pair (Carlesi et al. 2016). Halo assembly histories and the role of merger events have been investigated as well (Creasey et al. 2015), also finding that no recent major merger should have been occurred during the assembly of our MW (Scannapieco et al. 2015). It should be noted that although very limited in describing the details of baryonic physics, SAMs combined with DM simulations are an unavoidable tool to study the statistical properties of galaxy populations across a broad range of masses and redshifts (see for instance the recent *CATERPILLAR* project by Griffen et al. 2016a,b).

Despite a long series of investigations on how to implement mécanical and thermal feedback in different hydrodynamical schemes (Powell, Slyz & Devriendt 2011; Geen, Slyz & Devriendt 2013; Colin et al. 2016) and a continuous effort in performing code comparison projects (Scannapieco et al. 2012; Kim et al. 2014), we are still unable to consistently include feedback processes in models of galaxy formation and evolution. Feedback is often only partially implemented, even not understood in its basic physical principles, and depending on the problem at hand, a theoretical model would favour an accurate treatment of radiation transfer instead of a detailed gas dynamics, being their coupling feasible only under specific conditions. In other circumstances, for example when metal ions are used to trace metal enrichment, an alternative, detailed photoionization modelling is preferable (Graziani, Maselli & Ciardi 2013). Finally, some observable quantities can be simply more sensitive to photoionization than gas dynamics, depending on their physical time-scales.

In order to partially compensate these problems, the original version of the semi-analytic, data-constrained model of galaxy formation *GAMETE* ('GalaxyMERger Tree and Evolution'; Salvadori, Schneider & Ferrara 2007; Salvadori et al. 2010a) has been considerably extended (see de Bennassuti et al. 2014, hereafter DB14) and self-consistently coupled with the radiative transfer (RT) code *CRASH* ('Cosmological RAdiative transfer Scheme for Hydrodynamics'; Ciardi et al. 2001; Maselli, Ferrara & Ciardi 2003; Maselli, Ciardi & Kanekar 2009; Graziani, Maselli & Ciardi 2013) creating *GAMESH*, the first implementation of a full radiative feedback treatment in a semi-analytic model on top of a DM-only simulation. *GAMESH* is then capable to span the entire formation process of an MW-like halo across its cosmic evolution and to target many of the observable properties accessible in the local universe (see Graziani et al. 2015, hereafter LG15, for an introduction).

This paper is the first of a series where we will progressively exploit the many, new features of the *GAMESH* model and will apply it to a number of astrophysical problems, including the reionization history of the local Universe, the origin and spatial distribution of carbon-enhanced metal-poor stars in the MW halo and its satellites, the formation and coalescence environments of massive black hole binaries and the molecular and dust content of MW progenitors. Here we present how the *GAMESH* pipeline has been coupled to a new DM-only Galaxy formation simulation performed with the numerical scheme *GCD+* (Kawata & Gibson 2003), obtaining an isolated, MW-type halo in a cosmic cube of about 4 cMpc side length. The application of the SAM implemented in *GAMESH* on top of the new *N*-body simulation with an increased mass resolution allows us to perform an accurate analysis of the MW assembly history through accretion, mergers and dynamical interactions and to follow the redshift evolution of its baryonic properties, comparing them with observations at $z=0$. We investigate the chemical and star formation histories of MW progenitor galaxies and we critically compare them with observations of MW progenitor candidates in the redshift range $0.5 < z < 3$ (Papovich et al. 2015) and with observed scaling relations in $0 < z < 4$, including the galaxy main sequence (Schreiber et al. 2015), the stellar mass–metallicity relation, the fundamental metallicity and Fundamental Plane of metallicity relations (Mannucci et al. 2010; Hunt et al. 2012, 2016a).

The adoption of a first, simplified feedback scheme allows us to test the reliability of our galaxy formation model across redshift and to discuss the role that different galaxy progenitor populations play in the history of the MW and its Local Group (LG). Furthermore, it allows us to easily constrain the main free parameters of *GAMESH* in order to match the observed properties of the central MW-type halo avoiding complications introduced by the RT effects. Future investigations on the role of radiative and chemical feedback in galaxy formation, especially on low-mass galaxies, will adopt the RT described in LG15 and a new particle tagging scheme to simulate inhomogeneous metal enrichment. These features are essential to predict additional observable quantities, such as the MW metallicity distribution function and the luminosity function of its satellites.

The paper is organized as follows. In Section 2 we introduce the new DM simulation and we describe the halo catalogue, its assembly history and the properties of the resulting MW halo. A dedicated appendix also compares the new *N*-body simulation and the properties of the MW halo with four recent, independent simulations: *AQUARIUS* (Springel et al. 2008), *ELVIS* (Garrison-Kimmel et al. 2014), *CATERPILLAR* (Griffen et al. 2016b) and *APOSTLE* (Fattahi et al. 2016; Sawala et al. 2016); see Table A1 in Appendix A for more details.

<footnote n="1">¹ http://www.mpa.mpa-garching.mpg.de/millennium/</footnote>

<footnote n="2">² https://www.cosmosim.org/</footnote>
Section 3 describes the tuning of GAMESH on a set of observed MW properties. Section 4 focuses on analysing the properties of MW progenitors, while Section 5 analyses their evolution as hosted in minihaloes or Lyman $\alpha$ (Ly$\alpha$) cooling haloes. Section 6 finally summarizes the conclusions of the paper.

## 2 DM GALAXY FORMATION SIMULATION

In this section we describe the DM-only simulation performed to obtain an MW-size halo. We first describe the numerical scheme adopted in GCD+, the initial conditions of the simulation, the halo catalogue and its merger tree. The properties of the MW halo are finally described and the statistics of various halo populations in a surrounding volume of 4 cMpc side length. A careful comparison of the MW properties with similar haloes taken from independent simulations both in single and paired configurations can be found in Appendix A, where the halo properties are summarized in Table A2 for an easier comparison.

### 2.1 GCD+ and initial conditions

The N-body cosmological simulation of an MW-sized halo has been performed with GCD+ (Kawata & Gibson 2003; Kawata et al. 2013) with a β-version of periodic-boundary conditions and a TREEPM algorithm with parallel FFTW module. We used initial conditions created with MUSIC (Hahn & Abel 2011) and adopted a Planck 2013 cosmology ($\Omega_0 = 0.32, \Lambda_0 = 0.78, \Omega_b = 0.049$ and $h = 0.67$; Planck Collaboration XVI 2014) to simulate a volume of $83.5^3$ Mpc$^3$. In this volume we identified an MW-sized halo and we created the initial conditions for a zoom-in simulation. The final run consists of a total of 62 421 192 particles, 55 012 200 of which in the highest resolution region having particles with mass of $3.4 \times 10^5$ M$_\odot$. The virial mass of the resulting MW halo is $1.7 \times 10^{12}$ M$_\odot$.

To better resolve the early evolution of our universe we store the simulation outputs every 15 Myr from $z \sim 20$ down to $z = 10$, and every 100 Myr after this redshift. The total number of output snapshots is 155. This time resolution is high enough to follow the evolution of primordial stellar systems and to correctly account for gas recombinations through cosmic times, when the full pipeline mode including RT is adopted.

The final output of the simulation provides a list of collapsed halo objects and the projection on to grids of $512^3$ cells side$^{-1}$ of the DM distribution found in the 4 cMpc cosmic volume centred on the MW halo. From these grids the gas distribution in the cosmic web surrounding the MW is easily found by scaling the DM field with the value of the universal baryon fractions indicated by our cosmology. The resulting grid resolution in this domain is then $\sim 7.8$ ckpc.

### 2.2 Halo catalogue

We identify the populations of DM haloes at every snapshot by using a standard friend-of-friends (FoF) algorithm with a linking parameter of $b = 0.2$ and a threshold in the number of particles of 100. For each object we stored both virial properties (temperature $T_{vir}$, mass $M_{vir}$ and radius $R_{vir}$), as well the dynamical variables (position and velocity) of its centre of mass with respect to the central MW halo. Besides the list of DM haloes present in each snapshots, we also stored position and velocities of all particles resolving them. This information allows us to study the internal structure of the most massive haloes found in the LG (DM profiles, angular momentum, internal motion, overdensity structures, etc.), once they have a sufficient number of particles to reliably compute

these quantities. A detailed dynamical study of the MW satellite and subsatellite properties is still in preparation (see Mancini et al., in preparation). Before concluding this section, it is worth to note that the initial conditions of the simulation have not been selected to reproduce the structural and dynamical properties of the oLG, but rather to simulate a plausible MW-like halo and to focus on its evolution. As a result, the central 4 cMpc volume contains a total collapsed mass of $M_{DM} \sim 3 \times 10^{12}$ M$_\odot$ distributed in 2458 haloes. Among these, two haloes have DM mass $M \sim 10^{11}$ M$_\odot$ [M32, M33 or Large Magellanic Cloud (LMC)-type haloes, see table 1 in Guo et al. 2010 and references therein], 14 have $10^{10} \lesssim M < 10^{11}$ M$_\odot$ and 98 have $10^9 \lesssim M < 10^{10}$ M$_\odot$ (see Section 2.4 for their classification). Also note that the absence of a M31-like halo within 4 cMpc makes the total mass of our MW environment too low with respect to the oLG because $M_{M31} \gtrsim M_{MW}$ by recent estimates (see Ibata et al. 2007 and references therein). Other M31-sized haloes are found instead in the larger 8 cMpc region. Hereafter we will refer to this 4 cMpc cosmic region as the ‘LG’ of the present simulation. Also note that LG is also the maximum volume resolved exclusively by the high-resolution DM particles of GCD+ and also all its haloes are optimally resolved.

### 2.3 Merger trees and dynamical interactions

For each halo found by the FoF at redshift $z_i$, we have built its merger tree (MT) by iteratively searching all its particle IDs (pIDs) in the previous snapshots, back to the initial redshift $z_1$. A OpenMP$^3$ parallel searching technique, specifically tuned on the simulation data, has been developed to build-up the merger tree correlating the pIDs and halo IDs (hIDs) and to establish the ancestor/descendant relationships among hIDs found in $[z_i, z_{i-1}]$.⁴ It should be noted that once a pID at $z_{i-1}$ is not associated with any hID, it is associated with a reserved value we call ‘IGM ID’. We also verified that due to the recentring adopted for each snapshot to define the LG volume in the simulation data, few pIDs are sporadically not found in the LG volume because they are not geometrically captured. These pIDs are then classified as ‘missing’ and their associated hID marked as hID = -1, to exactly conserve the total mass. It is of primary importance to point out here that in our definition of merger tree the progenitors of a single halo do not necessarily collapse entirely into it, as generally assumed in semi-analytic tree models. In other words, in the merger trees of GAMESH descendant haloes do not conserve by design the total mass of their progenitors. The mass is instead conserved when contributions from progenitors are accounted for on a particle base.

While extremely demanding in term of computational processing, once done, this approach allows us to conserve the mass across dynamical interactions of DM haloes with their environments and to exactly follow all their dynamical processes regulating the accretion of DM haloes: mergers, tidal stripping and halo disruptions. All these events can then be classified and analysed and their baryonic counterpart accurately handled in the semi-analytic code. Baryonic properties (generally gas mass, metal mass and stars) can be then properly transferred throughout collapsed structures or returned to the baryonic intergalactic medium (IGM), by scaling with their relative dark mass contribution. While this is the approach used in the

<footnote>$^3$ www.openmp.org</footnote>

<footnote>$^4$ Note that once a particle belongs to many haloes, its multiplicity is also computed and stored to exactly account for the particle contribution in mass transfer processes across the merger tree steps.</footnote>
present paper, our scheme is sufficiently accurate to also associate halo baryons to the single DM particles in order to mimic an in-homogeneous spreading in the IGM of the LG by following DM particle dynamics. More details on the particle tagging scheme will be provided in future papers adopting it for specific applications. It should be noted that due to limitations introduced by the mass resolution of our DM simulation and the choice of our FoF parameters to identify haloes with 100 particles, our IGM could contain additional population of haloes with unresolved masses.

We classified the various dynamical interactions occurring during the mass assembly in the following categories.

(i) Halo growth by accretion: this event occurs when an isolated halo acquires particles only from the IGM, typically by mass accretion.

(ii) Halo growth by merger: when a halo at snapshot $z_{i+1}$ results in a contribution of two or many haloes at $z_i$, and possibly the IGM.

(iii) Halo stripping: when a halo loses part of its mass by tidal interactions with nearby haloes.

(iv) Halo destruction: when a halo found at $z_i$ loses its identity at $z_{i+1}$ because it is disrupted by tidal interactions and its particles are returned primarily in the IGM.

In the next section, we describe the assembly history of the most massive, MW-like, DM halo found at the centre of the LG cube at $z=0$.

## 2.4 The Milky Way DM halo assembly

Here we describe the assembly history of the MW halo defined above, in the context of its LG.

Fig. 1 shows the total mass of the MW merger tree (solid blue line) as a function of the lookback time ($t$). The dotted black line refers to the total collapsed mass found in the LG. As described above, our merger trees have been built particle-by-particle and we remind that the blue line shown here accounts, by design, for the total mass of the entire population of progenitor haloes providing particles collapsing on to the MW by the successive snapshot. As a consequence, the mass shown at certain $z_i$ does not necessarily transfer entirely to $z_{i+1}$ because of the complex series of dynamical interactions at play during the halo mass assembly.

To highlight the importance of having a complete merger tree we also show, as dashed red line, the MT resulting by following only the mass in most massive halo [$M_{MM}$, or major branch (MB)] at each $z_i$. It is immediately evident that by following the build-up history along the major branch, the discrepancy in mass becomes relevant at high redshift where a sensible fraction of the collapsed mass is distributed in a large number of MW progenitor haloes. The two MTs converge instead at $t = 4$ Gyr, (i.e. $z \sim 0.3$) where a large fraction of the final MW is already collapsed in $M_{MM}$ and all the progenitors contribute a large part of their mass. To understand the build-up time-scale, we computed the so-called ‘characteristic time’ ($t_a$) for the assembly of the MW halo, operationally defined as the redshift at which $M_{MM}(z) = M_{MW}(z=0)/2$ (see Mo et al. 2010) and it results in $t = 4.36$ Gyr, i.e. $z = 1.46$. Note that this estimate is compatible with the histories found in an independent set of DM simulations described in Behroozi, Wechsler & Conroy (2013, see in particular their fig. 19), and relative to haloes of similar mass. A further comparison with Table A2 shows that the MW halo has a formation redshift $z_{0.5}$ compatible with the ones of Cat-8 and Cat-12 (both isolated) while the extreme similarity with the formation redshift of Hamilton should be considered as a coincidence because this halo is found in a binary configuration and it is likely to have a very

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Build-up history of the central, MW-sized halo in the adopted N-body simulation. The total DM mass of all the MW progenitors ($M_{DM}$) is shown as function of the lookback time $t$ as solid blue line. The dashed red line shows the MW merger tree as obtained by following the major branch only. The total collapsed mass enclosed in the LG volume is shown as the dotted black line. For reference, the mass of similar MW-sized haloes taken from DM simulations or independent methods is also shown. Note that the scatter in these values is obtained from the scatter in the MW halo samples indicated in the various runs (see original papers for details). The inset panel illustrates the same MW history by plotting $M_{DM}(a)/M_{MW}(a=1)$ instead of the total collapsed mass.</figcaption></figure>

different accretion and dynamical history. Also note that the same time obtained from the blue line (i.e. accounting for all progenitors, see Navarro, Frenk & White 1996) results in $t = 2.66$ Gyr ($z \sim 2.45$), further indicating the importance of accounting for the large fraction of mass present as independent collapsed structures at high redshift.

Since in the literature the halo build-up histories are usually shown as $M_{DM}(a)/M_{MW}(a=1)$, in the bottom left-hand corner of this figure, we show the same history in these units for a more straightforward comparison with other simulations (see e.g. fig. 9 in Griffen et al. 2016b).

The resulting MW mass found at $z=0$ ($M_{MW} \sim 1.7 \times 10^{12} M_{\odot}$) is finally compared with the scatter in mass of MW-sized haloes<fnref n="5" /> found in recent simulations targeting the LG, both DM-only (the ELVIS simulation suite by Garrison-Kimmel et al. 2014 and the CATERPILLAR project by Griffen et al. 2016b) and the recent hydrodynamical APOSTLE simulations (Fattahi et al. 2016; Sawala et al. 2016); see Appendix A for more details. Finally note the additional agreement with the gold filled circle, showing estimates of Wang et al. (2015), obtained by using dynamical tracers.

To understand the MW growth within the global evolution of its LG, one can compare the blue and black solid lines of Fig. 1 and cross-check with the visual picture provided by Fig. 2, which shows the redshift evolution of the LG in a series of slice cuts intercepting

<footnote n="5">⁵ The mass ranges have been computed by extracting min/max mass values in tables of relative papers. Haloes in the scatter have been selected as 'reasonably close in mass' to our MW halo and just to suggest an indicative scatter introduced by different methods/simulations. More accurate details can be found both in Appendix A and in the original papers.</footnote>
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Slice cuts of the LG evolution at various redshifts. The panels show the DM density map obtaining by projecting the DM mass in each cell of the spatial grid. The total volume is 4 cMpc comoving mapped on a grid of 512 cells side⁻¹, for a spatial resolution of *r* ~ 7.8 ckpc.</figcaption></figure>

the central MW galaxy. Here the DM density map is shown as gradient from white (collapsed regions) to black (voids).<fnref n="6" />

It is immediately evident that while at high redshifts the mass of both MW and LG has a similar evolution,<fnref n="7" /> below *t* = 11 Gyr (*z* ~ 2) many structures not belonging to the MW merger tree, start collapsing in the entire volume or enter the domain from the larger scale.<fnref n="8" /> The evolution at high redshifts proceeds by assembling haloes along the diagonal filament created by the collapsing sheet. This is easily visible in the first slice cuts (top row, from left to right) where the time evolution of the main web filaments is shown. Below *z* ~ 3, the central halo dynamically dominates the LG region and continues to drag material entering from larger scales: around *z* ~ 2, an external filament not previously visible within 4 cMpc provides haloes to the central galaxy. This is a clear hint that galaxy formation is a multiscale process, assembling DM/baryonic mass created in different environments along the redshift (see also Section 4). At the final time (*z* = 0) the central halo shows a complex interplay with many filaments where a plethora of satellite galaxies are still collapsing towards the central attractor.<fnref n="9" />

As explained in Section 2.3, the accuracy of our MT allows us to disentangle the different growth processes (halo mergers or accretion from the IGM) and to describe their relative contribution. The result of this analysis is provided in Fig. 3 where we show the percentage of mass increase relative to the final MW mass ($\Delta M(t)/M_{\text{MW}}(z=0)$), as a function of the lookback time $t$ (solid red line). It is immediately evident that across cosmic times, the

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Differential contribution of DM mass, relative to the final MW mass, as a function of the lookback time $t$ (solid red line). The contributions from the IGM and collapsed structures are shown as dashed blue and green dotted lines, respectively.</figcaption></figure>

MW halo grows by means of a smooth and continuous assembly of matter spaced out by many violent accretion events, each of which provides a ~3 per cent contribution to the final mass [see e.g. the spikes around $t = 12.5$ Gyr ($z = 4.68$) and $t = 11.5$ Gyr ($z = 2.64$)]. As a further example note that a major event, increasing the mass of the most massive halo by about 5 per cent, is found around $t \sim 8$ Gyr (more precisely at $z = 0.95$), and this also corresponds to the last relevant major merger experienced by two Ly$\alpha$

<footnote n="6">⁶ Note that the equivalent gas number density is obtained by scaling the DM mass by the universal baryon fraction.</footnote>

<footnote n="7">⁷ This is mainly because all the haloes collapse first along a filament at the centre of the box.</footnote>

<footnote n="8">⁸ This is an 8 cMpc cube assumed to gravitationally constrain the LG domain and contains structures described by lower resolution particles.</footnote>

<footnote n="9">⁹ An animation can be found in the paper online resource files.</footnote>
cooling haloes ($T_{vir} > 2 \times 10^4$ K) found in the MW merger tree.<fnref n="10" /> Below $t = 3$ Gyr ($z \sim 0.2$) both the mergers and the accretion from the IGM phase become smoother and the mass growth progresses with steps contributing for less than 1 per cent to the final mass. Note that minor mergers are continuously found between the MW and small haloes or between $Ly\alpha$ cooling haloes and minihaloes ($T_{vir} \lesssim 2 \times 10^4$ K) orbiting the MW.

The relative contribution of accretion and mergers can be understood by comparing the dashed blue line (IGM accretion) and the dotted green line (halo mergers). While mass accretion from the surrounding IGM is dominating at all times, the green line shows an increasing number of halo mergers at low redshifts, with a substantial contribution around $t \in [2-4]$ Gyr (i.e. $z \in [0.15-0.35]$). It should be noted though, that while halo mergers contribute on average for less than 0.5 per cent to the final DM mass, their stellar, gas and metal contents contribute to shape the observed properties of the MW (see Section 3) and its surrounding satellite galaxies (see Section 4).

As pointed out in Wang et al. (2011), and also in agreement with the large samples of halo histories found in the Millennium-II simulation (Boylan-Kolchin et al. 2009), the growth of DM haloes can be dominated by mergers or by a smooth accretion of diffuse matter. Wang et al. (2011) consider six haloes of the *AQUARIUS* simulation (Springel et al. 2008) which target an MW-like halo mass at $z=0$ ($1-2 \times 10^{12}$ M$_\odot) and study their accretion mode and impact on the internal structure and age distributions of particles in the final halo. The authors claim that by averaging over the six haloes, smooth accretion can provide a relevant contribution to the final mass (roughly 30–40 per cent). We find that Aq-A-2 and Aq-C-2 experience a smooth growth history in their major branch, similar to our MW halo. When evaluated with a strict definition of mergers (1:3), the contribution of major mergers to the mass of Aq-A-2 and Aq-C-2 is below 0.1 per cent, in agreement with the estimates discussed above for MW.<fnref n="11" /> Interestingly, their mass and structural properties are also very similar (see details in Appendix A).

In this final paragraph we complement the information provided by the dynamics of the DM matter with the relative contribution of minihaloes and $Ly\alpha$ cooling haloes. In fact, the different impact of radiative and mechanical feedback on these populations strongly affects their evolution and leaves imprints on the observed properties of the MW and nearby dwarf galaxies, such as the metallicity distribution function of the most metal-poor stars in the Galactic halo (see Salvadori et al. 2010a; DB14; LG15; Salvadori et al. 2015; de Bennassuti et al. 2017, hereafter DB17).

Fig. 4 shows the redshift evolution of the number of minihaloes (thin lines) and $Ly\alpha$ cooling haloes (thick lines) along the merger tree of the MW (blue solid lines) and in the LG (black dotted lines). First note the plethora of minihaloes predicted by the DM simulation around the MW at $z=0$ ($N_{\text{Mini}} \sim 550$ within $2 \times R_{\text{vir}}$). In fact, during the last step of its mass assembly, less than 30 per cent of the entire minihalo population in the LG volume is embedded in the final MW halo; these haloes remain in the LG, providing a trace of the environmental conditions experienced by this volume along its redshift evolution. While the same considerations apply for the more massive population of $Ly\alpha$ cooling haloes, their number is about

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Number of haloes found in the LG (black dotted lines) and participating to the merger tree of the MW (blue solid lines) as a function of redshift. Minihaloes ($Ly\alpha$ cooling haloes) are shown with thin (thick) lines. The enclosed panel shows a zoom-in at high redshift.</figcaption></figure>

one order of magnitude lower ($\sim 3$ per cent), with very few objects ($N_{Ly\alpha} \sim 60$) still orbiting around the MW halo at $z=0$. The role of these two halo populations in shaping the baryonic properties of the central galaxy and its environment will be discussed in Section 5.

## 3 BARYONIC EVOLUTION OF THE MW

Here we follow the baryonic evolution of the galaxies associated with the DM haloes by running the semi-analytic part of the *GAMESH* pipeline. This means that radiative feedback is simulated by adopting a minimum mass of star-forming haloes and by assuming an instant reionization at $z_{\text{reio}} = 6$.<fnref n="12" /> A separate work will investigate the effects of radiative feedback on star-forming galaxies in the proper context of a local volume reionization simulation, performed by enabling the *CRASH* side of the pipeline (see LG15).

While more advanced versions of *GAMETE* introduce many improvements in the physical processes (Salvadori & Ferrara 2012; de Bennassuti et al. 2014, 2017; Salvadori et al. 2014) or adapt the original code to investigate the formation of quasars (Valiante et al. 2011, 2014, 2016; Pezzulli, Valiante & Schneider 2016), the pure semi-analytic scheme of *GAMESH* is based on the following simplifying assumptions.

(i) At each given time, stars are formed at a rate given by, star formation rate (SFR) = $\epsilon_*$ $M_{\text{gas}}/t_{\text{dyn}}$, where $\epsilon_*$ is the star formation efficiency, $M_{\text{gas}}$ is the total gas mass and $t_{\text{dyn}}$ is the dynamical time of the host halo.

(ii) Following Salvadori & Ferrara (2009, 2012), the star formation efficiency in minihaloes is assumed to be $\epsilon_{\text{MH}}/\epsilon_* = 2 \times [1 + (T_{\text{vir}}/(2 \times 10^4 \text{ K}))^{-3}]^{-1}$, as a result of the reduced efficiency of gas cooling (see also Valiante et al. 2016; DB17).

(iii) Stellar evolution is followed assuming the instantaneous recycling approximation (IRA). When the gas metallicity is

<footnote n="10">¹⁰ Note that this is not involving the most massive halo. In fact due to the peculiar history of our MW halo, its last major merger is found instead at $z > 5$.</footnote>

<footnote n="11">¹¹ Note that in this work we adopted a ratio (1:4) to identify a major merger event.</footnote>

<footnote n="12">¹² Alternative redshift modulated expressions can be found in Salvadori & Ferrara (2009, 2012), DB14, LG15 and DB17.</footnote>
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Redshift evolution of the total stellar (left-hand panel), gas (middle) and metal (right) masses. In each panel, the values predicted by GAMESH for the LG and for the merger tree of the MW are shown with black dotted and blue solid lines, respectively. The dashed green lines show the values computed by DB17 as averages over 50 independent Monte Carlo realizations of a semi-analytical merger tree for a $10^{12} M_{\odot}$ DM halo (the shaded areas show the corresponding $1\sigma$ deviation). Observations for the stellar and the gas masses at $z=0$ are taken from Flynn et al. (2006, green square), McMillan (2011, azure minus), Licquia & Newman (2015, violet downtriangle), Bovy & Rix (2013, yellow circle) and Kubryk, Prantzos & Athanassoula (2015, red uptriangle). The derivations of the mass of metals at $z=0$ are taken from Peeples et al. (2014, violet square) and DB14 (light green triangle). Note that in the middle and right-hand panels the blue solid and black dotted lines correspond to the mass of gas and metal used in star formation, while the green dashed lines with the shaded regions refer to the total gas and metal mass (see text).</figcaption></figure>

$Z < Z_{cr} = 10^{-4} Z_{\odot}$, Population III (Pop III) stars are formed with a constant mass of 200 $M_{\odot}$. Above the critical metallicity, Population II (Pop II) stars are formed with masses in the range [0.1–100] $M_{\odot}$, distributed according to a Larson initial mass function (IMF; Larson 1998) with a characteristic mass of $m_{\text{ch}} = 0.35 M_{\odot}$.

(iv) Chemical enrichment by supernovae and intermediate-mass stars is based on the same mass- and metallicity-dependent metal yields adopted in GAMETE and used in LG15.

(v) The mass outflow rate of supernova-driven winds is computed as $\dot{M}_{\text{gas,eje}} = 2 \epsilon_w v_{\text{circ}}^{-2} \dot{E}_{\text{SN}}$, where $\epsilon_w$ is the wind efficiency, $v_{\text{circ}}$ the host halo circular velocity and $\dot{E}_{\text{SN}}$ is the energy rate released by SN explosions, which depends on the SFR and on the stellar IMF (hence, a different value is adopted for Pop III and Pop II stars).

(vi) When $z \le z_{\text{reio}}$ star formation can only occur in galaxies with $T_{\text{vir}} > 2 \times 10^4$ K, to account for the effects of photoheating and photoevaporation (see LG15 for a thorough comparison between the instant reionization model and the model with a self-consistent reionization history computed by GAMESH).

We first implemented the simplest version of GAMETE in GAMESH because of many theoretical and practical (mostly numerical) reasons. First, previous runs with GAMETE on a semi-analytic merger tree (see DB17) and on top of a low-resolution N-body (see LG15) have shown that the simplest feedback implementation is sufficient to successfully calibrate the efficiency parameters of the model to fit the main integrated properties of the MW. This in turn, significantly reduces the range of possible values for our free parameters. Second, the introduction of a higher resolution simulation affects both the particle scheme of GAMESH and its RT side, so that the full pipeline becomes numerically demanding even on parallel facilities and then not suitable to make a blind parameter calibration. In future applications showing the new capabilities of the full RT scheme, we will add more observational constrains and will refine the calibration on the new set of observable quantities, also depending on the problem at hand.

The recalibration involves the two free parameters of GAMESH, namely the star formation efficiency in $Ly\alpha$ cooling haloes, $\epsilon_*$, and the efficiency of supernova-driven winds, $\epsilon_w$.

As discussed in LG15, we calibrate the free parameters of the model by requiring the SFR, the stellar and gas masses and the metallicity of the simulated MW galaxy at $z=0$ to match the observationally inferred values. For some of these quantities, such as the SFR or the total gas mass, the values inferred by different studies show up to one order of magnitude difference, as a result of the different tracers used in the observations or of the modelling strategy adopted to reconstruct the galaxy components (bulge, disc and halo). The interested reader can find in Kennicutt (1998), McKee & Ostriker (2007), Kennicutt & Evans (2012) and Bland-Hawthorn & Gerhard (2016) a large collection of critically revised estimates and galaxy modelling techniques.

Besides the different methodologies, in the last years the total stellar mass ($M_*$) inferred for the MW has largely converged to a value of $M_* = [3-7] \times 10^{10} M_{\odot}$, as proven by a series of independent estimates (Flynn et al. 2006; McMillan 2011; Bovy & Rix 2013; Kubryk et al. 2015; Licquia & Newman 2015). These are shown by coloured points in the left-hand panel of Fig. 5. In the same panel, we show the redshift evolution of $M_*$, as predicted by the more advanced SAM described in DB17 and ran on top of a semi-analytic merger tree (solid green line with the shaded region) and the stellar mass assembly predicted by GAMESH in the LG and for the merger tree of the MW, when $\epsilon_* = 0.09$ and $\epsilon_w = 0.0016$. The model is in good agreement with the observations, with a final value of $M_* \sim 4.6 \times 10^{10} M_{\odot}$ for the MW halo candidate. It predicts a total stellar mass of $\sim 6 \times 10^{10} M_{\odot}$ in the LG and its redshift evolution results consistent with that predicted by DB17 and its statistical scatter. A word of caution is also necessary here when interpreting the evolution in redshift of the baryons accounted for in the GAMESH merger tree and shown in this figure. While GAMESH transfers the baryons from progenitors to descendants exactly scaling by the DM particle contribution, along the merger tree lines shown here the masses do not conserve in redshift, as commented in the DM evolution session. At fixed redshift the mass shown in the merger tree line is the total baryonic mass of the progenitors and not their contribution in mass to the descendants. It is then an estimate of the maximum potential mass available from halo progenitors and the accretion from IGM. Also remember that the mass shown in
the LG does not conserve across redshift because of the continuous exchange of systems with the larger scale.

Below $z \sim 4$ the models show a different evolution, and at $z=0$ the DB17 predicts a mass $M_* \sim 7 \times 10^{10} M_{\odot}$ with a local star formation efficiency $\epsilon_* = 0.8$, i.e. one order of magnitude higher than the value required by *GAMESH*.<fnref n="13" /> The reason for this difference can be ascribed to (i) the different mass of the final MW DM halo, which in DB17 is assumed to be $M_{MW}(z=0) = 10^{12} M_{\odot}$, a factor of 1.7 smaller than the value assumed by the *N*-body simulation adopted here; (ii) the different DM evolution of the MW halo predicted by the *N*-body simulation, with continuous mass exchanges between haloes entering the MW merger tree and haloes of the LG. Conversely, the semi-analytic merger trees (based on the extended Press–Schechter formalism) are, by construction, mass conserving, so that any stellar population formed at $z > 0$ along the merger tree will inevitably end up in the MW by $z=0$; (iii) the presence of the IRA that accelerates stellar evolution and underestimates the mass of active stars present at each given time.

In addition, at $z \le 4$ the *GAMESH* MW halo progenitors have systematically one order of magnitude higher total gas mass than the ones predicted by the semi-analytic merger tree adopted in DB17.

This can be seen in the middle panel of Fig. 5, where we show the redshift evolution of the gas mass. Here the black dotted and blue solid lines represent the mass of gas used in star formation, i.e. $\epsilon_* M_{gas}$, while the green dashed line represents to the total $M_{gas}$ predicted by DB17. As illustrated in Section 2.4 and in Fig. 2, the early assembly of the MW halo in the *N*-body simulation is dominated by mass accretion and mergers of nearby haloes, and the evolution is similar to the one predicted by DB17 using the semi-analytic merger trees based on the extended Press–Schechter (EPS) formalism (Press & Schechter 1974). When $z < 4$, the simulated MW halo grows by many episodes of violent accretion and many minor mergers with haloes entering the LG from the larger scales. All these effects cannot be accounted for by the semi-analytic merger trees.

Despite the intrinsic differences found in their assembly histories, both models predict a final gas mass in the MW in agreement with the observed values. Here the comparison among models and with observations should be taken with caution. In fact, the total gas mass in the MW includes the cold and warm components (molecular and atomic phases mainly in the disc) and a hot halo (coronal) component, as exhaustively detailed in Bland-Hawthorn & Gerhard (2016). The observations reported in the middle panel of Fig. 5 refer to the total interstellar medium (ISM) mass in the disc as inferred from dynamical measurements (Bovy & Rix 2013, yellow circle), and by averaging the values of the atomic and molecular gas masses obtained by different observational studies (Kubryk et al. 2015, red triangle). We note that by adopting the most likely mass range for the Galactic corona ($2.5 \pm 1) \times 10^{10} M_{\odot}$ (Bland-Hawthorn & Gerhard 2016), the total baryonic mass is estimated to be in the range $[7-11] \times 10^{10} M_{\odot}$. If we account for the total amount of gas enclosed in the MW halo, we find $M_{gas} \sim 1.3 \times 10^{11} M_{\odot}$.

Finally, in the right-hand panel of Fig. 5 we show the evolution of the metal mass in the ISM. Here we do not follow separately the evolution of dust, hence all the lines show the total mass in heavy elements (gas-phase metals and dust). The models show behaviours

that reflect their corresponding stellar mass assembly histories, and DB17 predicts a higher metal content at $z=0$, consistent with its higher $M_*(z=0)$. Both models are in agreement with the violet square at $z=0$, based on a detailed inventory of metal mass components in present-day $L_*$ galaxies (Peeples et al. 2014). To be consistent with the observations reported in the other panels of Fig. 5, we have computed the mass of metals and dust in the ISM from the fitting functions of Peebles et al. (2014), using a stellar mass in the range $[3-7] \times 10^{10} M_{\odot}$ and we find $M_Z = [0.95-4.7] \times 10^8 M_{\odot}$.

In addition, Peebles et al. (2014) provide an estimate of the circumgalactic medium (CGM) metal mass as probed by low- and high-ionization (O VI) species based on Cosmic Origins Spectrograph (COS)-Halos data and of the mass of dust in the CGM based on the reddening of background quasars (Ménard et al. 2010). Taking the values from their table 5 and adding the metal mass in the hot X-ray emitting CGM gas,<fnref n="14" /> we infer a total mass of heavy elements in the CGM of $1.7 \times 10^8 M_{\odot}$ with minimum (maximum) values of $0.9(3.7) \times 10^8 M_{\odot}$. Hence, the total estimated mass of metals is found to be in the range $[1.85-8.4] \times 10^8 M_{\odot}$. Our model predicts a total mass of metals in the final MW halo to be $\sim 10^9 M_{\odot}$. However, at the MW mass scale, the observations probe the metal content within 150 kpc (Peeples et al. 2014). If we assume the metals to follow the DM halo radial profile at these large radii, we find $M_Z(r < 150 \text{ kpc})$ predicted by the simulation to be $6.9 \times 10^8 M_{\odot}$, in agreement with the observed value.

Hence, we conclude that having selected the two free parameters to be $\epsilon_* = 0.09$ and $\epsilon_W = 0.0016$ the stellar mass and the mass of gas and metals in the ISM predicted by *GAMESH* for the MW-like halo at $z=0$ are consistent with the observations. A multiphase treatment of the gas and dust evolution in *GAMESH*, similar to what presented in DB14 and DB17, and the relaxation of the IRA will be implemented in the future to address specific problems and are not expected to sensitively change the previous calibration.

In Fig. 6 we show the total SFR predicted by our model as function of the lookback time. Data points indicate observationally inferred estimates found in the literature, from the oldest estimates by Smith, Biermann & Mezger (1978, SFR ~5 $M_{\odot}$ yr$^{-1}$, dark green cross), McKee & Williams (1997) and Diehl et al. (2006, SFR ~4 $M_{\odot}$ yr$^{-1}$, violet square) to the newest, generally lower values around SFR ~2 $M_{\odot}$ yr$^{-1}$ by Robitaille & Whitney (2010, orange minus), Licquia & Newman (2015, light-violet circle) and Kubryk et al. (2015, red triangle).

As exhaustively discussed in Chomiuk & Povich (2011), the number of assumptions needed to derive the total SFR of the current MW (e.g. in its structure, its stellar sample and stellar IMF) is so large that the resulting scatter can span one order of magnitude (see the azure diamond with the largest error bars). For the MW-like halo at $z=0$, *GAMESH* finds a SFR ~4.7 $M_{\odot}$ yr$^{-1}$, a factor of 2 higher than recent estimates but still compatible with the data scatter. In the enclosed panel, we show the same SFR as a function of redshift $z$, to better visualize the evolution at high redshift.

While the global trend of the SFR predicted by *GAMESH* agrees with the one found in DB17, and the two SFRs show similar peak values (~15–20 $M_{\odot}$ yr$^{-1}$), they peak at different redshifts. The progressive, quasi-parallel decline of the total SFR results in final SFRs differing by ~1.5 $M_{\odot}$ yr$^{-1}$. As argued for the gas and stellar mass behaviours, we ascribe these discrepancies to intrinsic differences in the merger tree definitions discussed above, MW mass

<footnote n="13">¹³ Note that in both models the global star formation efficiency is defined as $M_*/M_{gas} = \epsilon_*(\Delta t/t_{dyn})$, where $t_{dyn}$ is the redshift-dependent halo dynamical time and $\Delta t$ is the time-scale assumed for star formation.</footnote>

<footnote n="14">¹⁴ We compute this quantity from equation (24) in Peebles et al. (2014) assuming a stellar mass of $[3-7] \times 10^{10} M_{\odot}$.</footnote>
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Total SFR in the MW merger tree (solid blue line) and in the LG (dotted black line) as function of the lookback time. The average SFR found in DB17 is shown with dashed green line, with the shaded region showing the $1\sigma$ dispersion. In the enclosed panel the same quantities are shown as function of $z$. Data points with error bars, when available, are taken from the literature (see text for details).</figcaption></figure>

assembly history, particularly at $z < 4$, and to the IRA, which naturally introduces an acceleration in the evolution.

A comparison of *GAMESH* and its predictions with independent SAMs can be found in Appendix B.

## 4 PROPERTIES OF MILKY WAY PROGENITORS

So far, we have investigated the global properties of the simulated haloes, in the MW merger tree and in the LG. In this section, we discuss the SFR, the mass in stars, gas and metals predicted for MW progenitor systems at $0 < z < 4$ by the *GAMESH* simulations and compare these with observations.

Recent studies have started to investigate the redshift evolution of progenitors of MW-like galaxies at $z=0$, selecting candidates from very deep near-infrared (IR) surveys on the basis of their constant comoving density (van Dokkum et al. 2013), of their evolution

on the galaxy star-forming main sequence (Patel et al. 2013) or of multi-epoch abundance matching techniques (Papovich et al. 2015).

Using a combined data set based on the FourStar Galaxy Evolution (ZFOURGE) survey, Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS), *Hubble Space Telescope (HST)*, *Spitzer* and *Herschel*, Papovich et al. (2015) derived photometric redshifts and stellar masses for MW progenitors and discuss their evolution with redshift. To compare with their analysis, we have extracted MW progenitors from the *GAMESH* simulation adopting a similar selection procedure. We first identify all the simulated systems in the same redshift bins of the observations, and then we select those with a stellar mass that falls within ±0.25 dex of the central stellar mass adopted by Papovich et al. (2015, see entries 1 and 2 in their table 1). The number of selected progenitors ranges between 41 (in the lowest redshift bin, 0.2 < z < 0.7) and 6 (in the highest redshift bin, 2.2 < z < 2.8).

Fig. 7 shows the resulting evolution of the average stellar mass (left-hand panel), SFR (middle panel) and gas fraction (right-hand panel), defined as $f_{\text{gas}} = M_{*}/(M_{\text{gas}} + M_{*})$. The blue solid lines show the model predictions, with the shaded region representing the $1\sigma$ scatter, and the red points are the Papovich et al. (2015) data. To increase the statistics of MW progenitors, particularly in the higher $z$ bins, we also show the model predictions when the mass selection is done within ±1 dex of the central mass adopted by Papovich et al. (2015, grey dotted line with shaded region). Our simulation suggests that MW progenitors follow a stellar mass trend that is in good agreement with the observations, particularly if the mass selection includes a larger number of MW progenitor systems at $z > 1.5$. In agreement with previous studies, we find that more than 90 per cent of the MW mass has been built since $z \sim 2.5$. However, the SFR and the gas fraction of the simulated galaxies have a shallower evolution in the 3 Gyr period between $z = 2.5$ and 1 than found by Papovich et al. (2015). In particular, the peak SFR of $\sim 10 M_{\odot} \text{yr}^{-1}$ at $z \sim 1-2$ of the most massive MW progenitors is smaller than the value reported by Papovich et al. (2015), and in closer agreement with the evolution found by van Dokkum et al. (2013, see the red dashed lines). We find that the MW mass build-up can be fully explained by the SFRs of its progenitor systems, and does not require significant merging (van Dokkum et al. 2013). If star formation dominated the formation of the MW galaxy, then its growth must heavily depend on the evolution of cold gas and gas-accretion histories. This is consistent with the results presented in Section 2. In addition, by inverting the Kennicutt–Schmidt law, Papovich et al. (2015) show

<figure><img src="image_8.png" /><figcaption>**Figure 7.** Comparison between the stellar mass (left-hand panel), SFR (middle panel) and gas fraction (right-hand panel) of MW progenitors in the *GAMESH* simulation and observational data from Papovich et al. (2015, red points). The blue solid lines show the average values among MW progenitors selected in each redshift bin following the same procedure adopted by Papovich et al. (2015). The grey dotted lines show the results when the minimum mass adopted to select MW progenitors in each redshift bin is decreased by 1 dex (see text). The shaded regions represent the $1\sigma$ scatter around the mean.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 8.** The SFR as a function of stellar mass of all MW progenitors in different redshift bins, as indicated in the legend. In each panel, the points represent the simulated systems and the dashed line shows the analytic fit to the galaxy main sequence at the central redshift of the bin, taken from Schreiber et al. (2015). The dotted lines are a factor of 2 above/below the fit.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** The mass-metallicity relation at different redshifts (see the legends). The points show the simulated MW progenitors and the solid lines represent the fit to the observed relations reported by Maiolino et al. (2008, at $z < 3$) and Mannucci et al. (2009, at $3 < z < 4$). The dashed lines are the fit obtained by Hunt et al. (2016a) with the shaded region showing the ±0.15 dex scatter. The tight relations followed by some of the points, particularly in the top right-hand panel, identify evolutionary tracks of galaxies in the redshift range encompassed by the observed samples.</figcaption></figure>

that the effective size and SFRs imply that the baryonic cold-gas fractions drop as galaxies evolve from high redshift to $z \sim 0$ (see the red data points in the right-hand panel of Fig. 7). The predicted $f_{\text{gas}}$ of the simulated sample show instead a rather flat trend and, independently of the adopted selection criteria, the average SFR and gas fraction are larger than inferred by the observations below

$z \sim 1$. The above-mentioned discrepancies in the evolution of the average SFR and gas fraction might be induced by the assumed IRA, which affects the efficiency of mechanical and chemical feedback acting on the evolution of individual galaxies, and by the lack of RT effects. All these different feedback processes, indeed, can strongly affect the evolution of the MW progenitors (e.g. Salvadori, Dayal & Ferrara 2010b; Salvadori & Ferrara 2012; LG15).

We have also checked the position of the simulated MW progenitors relative to the galaxy main sequence of star formation. In Fig. 8 we show the results using the same redshift bins adopted in Fig. 7, but without making any selection on the stellar mass. In each panel, the points represent all the simulated systems, while the dashed line is the analytic fit to the observations, taken from Schreiber et al. (2015, see their equation 9) and computed at the central redshift of each bin. There is a large scatter in the SFR of the smallest MW progenitors and most of the systems with $M_* < 10^8 \text{ M}_\odot$ show SFRs that can vary by almost one order of magnitude. While the galaxy main sequence cannot be constrained by observations in this regime, an increasing scatter towards low stellar masses has already been found in hydrodynamical simulations as a result of the rising importance of stellar feedback (Hopkins et al. 2014). Yet, the more massive among the MW progenitors at each redshift lie within a factor of 2 of the galaxy main sequence (the region within the two dotted lines) all the way from $z \sim 2.5$ to $\sim 0$.

Finally, we compare the gas metallicity of the simulated MW progenitors with the observed mass-metallicity relation (MZR) at different redshifts and with two (redshift independent) combinations of stellar mass, SFR and metallicity known as the fundamental metallicity relation (FMR; Mannucci et al. 2010) and Fundamental Plane of metallicity (FPZ; Hunt et al. 2012, 2016a).

The results are shown in Figs 9 and 10, respectively. The interstellar oxygen abundance has been computed assuming a solar oxygen-to-metal mass ratio of 0.00674/0.0153 = 0.44 (Caffau et al. 2011), so that the solar metallicity corresponds to $12 + \log(O/H) = 8.759$. The points represent all the simulated MW progenitors in the same redshift bins of Maiolino et al. (2008) and Mannucci et al. (2009), without any additional selection on the stellar mass or SFR. Instead, the solid lines show the fit to the data and are drawn only in the mass range probed by the observations. It is clear that most of the simulated systems have stellar masses that are outside this range, except for the few most massive MW progenitors with $M_* > 10^9 \text{ M}_\odot$. Similarly to Fig. 8, the scatter in the MZR increases with decreasing $M_*$, a result that appears to be consistent with deep spectroscopic observations that probe galaxies down to $M_* \sim 3 \times 10^7 \text{ M}_\odot$ at $0.5 \le z \le 0.7$ (Guo et al. 2016). At $3 < z < 4$, the most massive systems have metallicity slightly higher than those implied by the Mannucci et al. (2009) fit. However, at $z < 3$ the simulated systems fall systematically below the fits by Maiolino et al. (2008). A better agreement is found with the fit to the MZR proposed by Hunt et al. (2016a, and computed using their equation 2 at the average redshift of each bin), shown as the dark grey dashed line, with the shaded region representing a dispersion of ±0.15 dex.

A similar result is found in Fig. 10, where we show the position of the simulated MW progenitors relative to the FMR (left-hand panel) and FPZ (right-hand panel). Here we have reported systems selected in the same redshift bins as in Fig. 9. As usual, smaller MW progenitors that populate the lower right-hand side of the panels show a large scatter at all redshifts. The most massive MW progenitors align along the FMR but with a −0.5 dex metallicity offset. Conversely, their metallicity is within the scatter of the FPZ.

We conclude that while the simulated systems may be slightly too metal-poor at high stellar masses and too metal-rich at lower
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Distribution of the MW progenitors relative to the fundamental metallicity relation by Mannucci et al. (2010, left-hand panel) and to the Fundamental Plane of metallicity by Hunt et al. (2016a, right-hand panel). All the simulated systems at $\langle z \rangle = 0.07, 0.7, 2.2$ and at $3 < z < 4$ (using the same colour coding adopted in Fig. 9) are shown with data points. Dashed lines with shaded regions show the observed fits (see text).</figcaption></figure>

stellar masses, the discrepancy with the MZR evolution by Maiolino et al. (2008) and with the FMR by Mannucci et al. (2010) may be at least partly due to the different metallicity calibrations used by these authors, which may overestimate the observed metallicity by 0.3 dex (Hunt et al. 2016a, see in particular their section 4.4). (More details about discrepancies in metallicity calibrations can be found in Kewley & Ellison 2008.)

More interestingly, the distribution of the most massive MW progenitors is consistent with the FPZ and aligned with the FMR. Since these redshift-independent scaling relations between metallicity, stellar mass and SFR are believed to originate from the interplay between gas accretion, star formation and SN-driven outflows (see, among others, Dayal, Ferrara & Dunlop 2013; Hunt et al. 2016b), we conclude that the description of these physical processes in the GAMESH simulation leads to results consistent with observations at $0 < z < 4$ even in a simulated, biased region of the current Universe.

## 5 EVOLUTION OF GALAXIES HOSTED IN MINIHALOES AND LY$\alpha$ COOLING HALOES

In this section, we study the evolution of galaxy populations by assessing their relative contribution to the total quantities discussed above. By following the standard classification of DM haloes introduced in Section 2 (see Fig. 4), we discuss the properties of the galaxies hosted in minihaloes and Ly$\alpha$ cooling haloes. We first follow their SFR, $M_*$, $M_Z$ in time and then we compare their final distribution at $z=0$.

In Fig. 11 we show the redshift evolution of their SFR, both in the MW merger tree and in the LG. In the left-hand panel, the total SFR of minihaloes and Ly$\alpha$ cooling haloes found in the MW merger tree is shown by solid thin and thick blue lines, respectively. At all but the highest redshifts ($z>17$), the MW SFR is dominated by galaxies hosted in Ly$\alpha$ cooling haloes.<fnref n="15" /> Despite the number of minihaloes

is largely dominant, as shown in Fig. 4, their small contribution to the total SFR at all redshifts is due to (i) their low star formation efficiency (Salvadori & Ferrara 2009, 2012; LG15) compared to Ly$\alpha$ cooling haloes, where the SFR is proportional to the halo gas content and controlled by the accretion history (see Fig. 3); (ii) their low gravitational potential, which implies an intrinsic inefficiency in accreting gas from the IGM.

Among Ly$\alpha$ cooling haloes, a special role is played by the most massive halo that drives the major branch of the MW. The SFR along the MB is shown by the red dashed line. The comparison with the solid blue line shows that the MB dominates the SFR in the MW only at $z < 1$. This does not come as a surprise, given that the MB contains half of the mass of the final MW halo at $z \lesssim 1.46$ and that at higher $z$ the MW mass assembly is dominated by a multiplicity of Ly$\alpha$ cooling systems, which also dominate its SFR. This is clearly visible in the inset of the left-hand panel of Fig. 11, where we show the SFR in the MB as a function of the lookback time (dashed red line) and the SFR of individual Ly$\alpha$ cooling haloes that are part of the MW merger tree but not collapsed in the MB (blue points). The progressive disappearance of these points is a consequence of their accretion on to the MB. We note that the flat evolution of the SFR in the MB in the past 8 Gyr is in line with findings of independent models (De Lucia et al. 2014).

In the right-hand panel of Fig. 11, we investigate the relative contribution of Pop II and Pop III stars to the SF history in the merger tree of the MW and in the LG. In the very early evolution ($z>16$), the SFR is dominated by Pop III stars in both the LG and along the MW merger tree (represented by the cyan and blue shaded areas, respectively). Due to the effect of metal enrichment, their formation is mostly confined to the first star-forming minihaloes and the sharp drop in their SFR is driven by the increase of the average metallicity of the IGM above $Z_{crit} = 10^{-4} Z_{\odot}$ at $z < 16$. We note that the Pop III SFR is larger in the LG than along the MW merger tree, as a result of the larger statistics of minihaloes. This suggests that traces of early Pop III star formation are not confined to the MW and its satellites but may be found in external galaxies of the LG, although their detectability seems to remain beyond the capabilities of the new generation of telescopes, such

<footnote n="15">¹⁵ The curves that show the contribution of minihaloes are interrupted at $z=6$ because we assume that SF is suppressed in minihaloes when $z \le z_{reio}$ (see Section 3).</footnote>
<figure><img src="image_12.png" /><figcaption>**Figure 11.** Redshift evolution of the total SFR of galaxies hosted in minihaloes and Lyα cooling haloes along the merger tree of the MW and in the LG. Left-hand panel: SFR of minihaloes and Lyα cooling haloes belonging to the MW merger tree (solid thin and thick blue lines, respectively). The dashed red line shows the SFR along the major branch of the MW (the most massive halo at each redshift). The panel inset shows the SFR along the major branch (red dashed line) as function of the lookback time $t$ (Gyr) and the SFR of all the Lyα cooling haloes as blue points. Right-hand panel: SFR of Pop II stars along the MW merger tree and hosted in minihaloes (solid thin red line) and in Lyα cooling haloes (solid thick red line). The dotted black line indicates the Pop II SFR history in all the minihaloes of the LG. Pop III SFRs along the MW merger tree and in the LG are indicated by shaded areas (blue and cyan, respectively). For comparison, we also show the SFR along the MW MB (red dashed line). The panel inset shows the SFR versus stellar mass of all galaxies hosted by Lyα cooling haloes in the LG at $z=0$.</figcaption></figure>

as the James Webb Space Telescope (JWST).<fnref n="16" /> Finally note how the sharp transition between Pop III and Pop II stars is certainly affected by the absence in the current study of inhomogeneous RT and metal mixing in the LG, as predicted by independent models (Scannapieco, Schneider & Ferrara 2003; Salvadori et al. 2014; Sarmento, Scannapieco & Pan 2017). This point will be deeply investigated in the companion work of this paper including full RT and particle tagging (Graziani et al., in preparation).

The remaining redshift evolution is dominated by Pop II stars formed in minihaloes and Lyα cooling haloes (solid thin and thick red lines, respectively). We also show the contribution of Pop II star formation in minihaloes belonging to the LG (black dotted line). It is interesting to note that this component is comparable to the total SFR along the merger tree of the MW and it is larger than the SFR on the MB down to $z \approx z_{reio}$, below which it is suppressed by reionization. Although our simplified description of radiative feedback may overestimate star formation in small objects prior to reionization, this comparison shows that they provide an important source of ionizing photons within the LG (Salvadori et al. 2014; Graziani et al., in preparation).

The star-forming systems found in the LG at $z=0$ are shown in the figure inset, where their SFR is plotted as function of the stellar mass.<fnref n="17" /> Their distribution follows the galaxy main sequence presented in Fig. 8. The three most massive star-forming galaxies, with $M_* > 10^9 M_\odot$, have $0.1 < \text{SFR} < 1 M_\odot \text{ yr}^{-1}$, similar to the values typically found in big objects observed in the LG.

Similarly, galaxies with lower stellar mass, $10^8 < M_*/M_\odot < 10^9$ have $10^{-2} < \text{SFR} < 10^{-1} M_\odot \text{ yr}^{-1}$, as observed in the LMC and Small Magellanic Cloud (SMC) that have $M_* \approx 10^9 M_\odot$, $0.25 < \text{SFR} < 0.63 M_\odot \text{ yr}^{-1}$ (LMC) and $M \approx 3 \times 10^8 M_\odot$, $0.016 < \text{SFR} < 0.039 M_\odot \text{ yr}^{-1}$ (SMC; Skibba et al. 2012). Our simulation also finds a third population of galaxies with $M_* < 10^8 M_\odot$ hosted in small Lyα cooling haloes with SFR $< 10^{-3} M_\odot \text{ yr}^{-1}$, within the detection limits of recent dwarf galaxy surveys, such as the Dwarf Galaxy Survey (DGS; Madden et al. 2013, 2014).

In the two panels of Fig. 12 we show the redshift evolution of the stellar (top panel) and metal (bottom panel) mass, separating the contribution of MW progenitors hosted in minihaloes and Lyα cooling haloes (solid thin and thick blue lines, respectively), and of all the minihaloes found in the LG (dotted black lines). As expected, progenitors in Lyα cooling haloes dominate the evolution and the total stellar and metal mass in minihaloes is orders of magnitude smaller, even when all the minihaloes in the LG are considered (with the exception of the metal mass at $z > 16$, due to efficient metal enrichment by Pop III stars). The evolution of the stellar mass of both Lyα cooling and minihaloes is in good agreement with the results of Salvadori et al. (2014) that match the observed properties of MW dwarf galaxies using a more sophisticated version of the SAM in which the IRA approximation is relaxed and the inhomogeneous metal enrichment and reionization of the MW environment are accounted for by using analytic prescriptions.

At $z < 6$, the stellar mass in minihaloes decreases as a result of the combined effect of radiative feedback following reionization and assembly of minihaloes into bigger structures (see also Fig. 4).

Below $z < 2$, the oscillating behaviour in the thin blue lines is due to continuous mass exchange between small progenitors (both mini and small Lyα cooling haloes) orbiting around the MB, before the

<footnote n="16">¹⁶ Due to the intrinsic differences in the dynamical configuration of haloes in our LG with respect to the observed one, this information should be valued only in a statistical sense.</footnote>

<footnote n="17">¹⁷ Note that the MW galaxy is excluded.</footnote>
<figure><img src="image_13.png" /><figcaption>**Figure 12.** Contribution of galaxies hosted in minihaloes and Ly$\alpha$ cooling haloes to the total mass in stars (top panel) and metals (bottom panel) as a function of redshift. MW progenitors hosted by minihaloes and Ly$\alpha$ cooling are shown by the solid thin and thick blue lines. The dotted black lines show the systems hosted in minihaloes in the LG.</figcaption></figure>

final merging, which causes the sharp drop at $z = 0$.<fnref n="18" /> A series of mass exchange events through dynamical interactions or destructive mergers involving small Ly$\alpha$ cooling and minihaloes is at the origin of the final, spiky mass increase in the dotted black lines, referring to small dwarf galaxies surviving in the LG at $z = 0$. By checking the merger trees of these small objects surviving in the local Universe, we verified that the increase in metallicity of the minihalo population is mainly ascribable to a series of contamination events from highly enriched baryon masses tidally dragged by minihaloes orbiting Ly$\alpha$ cooling haloes. An increasing number of small Ly$\alpha$ cooling haloes transitioning the minihaloes population after mass loss by tidal interactions is also found in the lowest redshift evolution. A deeper understanding of this interesting interplay between systems dominated by interacting dynamics and sharing tracers of chemical evolution, certainly requires the adoption of a full particle scheme and it is the deferred to a future work.

# 6 CONCLUSIONS

In this paper we explored the properties of MW progenitors by simulating the Galaxy formation process with the GAMESH pipeline. To guarantee a good statistics of halo progenitors, GAMESH adopted a new DM simulation providing a well mass resolved cosmological box containing a central MW-like halo and a significant number of smaller objects and MW satellites. This new simulation allows us to draw the following main conclusions.

(i) The new DM simulation has the adequate mass resolution to guarantee a reliable MW-type DM halo whose mass and dynamical properties are in agreement with a number of independent simulations and observational estimates. The LG surrounding it shows a plethora of intermediate-mass Ly$\alpha$ cooling haloes and a vast number of minihaloes, useful to study both the global accretion process and the effects of mergers and tidal stripping. The time resolution of the new simulation has proven to be adequate to reproduce the major events of the accretion history of the MW halo, in agreement with similar trends found by independent simulations.

(ii) Once processed by GAMESH and after tuning the star formation and galactic wind efficiency, the stellar gas and metal mass at $z = 0$ are consistent with MW observations. We predict a final SFR for the resulting MW system a factor of 2 higher than recent simulations but still in agreement with the many observational uncertainties.

(iii) A particle-by-particle reconstruction of the MW merger tree allows us to follow the redshift evolution of MW progenitor galaxies and to predict their baryonic properties. The simulated progenitor galaxies follow a stellar mass trend in good agreement with observations targeting 'plausible' MW progenitors (van Dokkum et al. 2013; Papovich et al. 2015). Our simulation suggests that more than 90 per cent of the MW mass has been built since $z \sim 2.5$. However, the SFR and the gas fraction of the simulated galaxies have a shallower evolution between $z = 2.5$ and 1 than found by Papovich et al. (2015). While the MW mass build-up can be fully explained by the SFRs of its progenitor systems, and does not require significant merging, we will re-evaluate all the discrepancies found with observations in future model implementations relaxing the IRA approximation and accounting for detailed radiative feedback.

(iv) The most massive among the simulated MW progenitors lie within a factor of 2 of the galaxy main sequence all the way from $z \sim 2.5$ to $z \sim 0$. The predicted SFRs show an increasing scatter towards low stellar mass systems due to the rising importance of feedback effects. Similar results are found when comparing the distribution of the simulated galaxies with the observed MZR, FMR and FPZ at $0 < z < 4$ (Mannucci et al. 2010; Hunt et al. 2012, 2016a). Since these scaling relations are believed to originate from the interplay between gas accretion, star formation and supernova-driven outflows, we conclude that the description of these physical processes obtained by GAMESH leads to results consistent with observations.

(v) At all but the highest redshifts, the SFR of the MW is dominated by a multiplicity of galaxies hosted in Ly$\alpha$ cooling haloes, hosting Pop II stars. These systems are progressively accreted by the major branch of the MW merger tree, which provides the dominant contribution to the SFR at $z < 1$. The cumulative contribution of star-forming minihaloes in the LG is comparable to the SFR along the MW merger tree at $z > 6$, indicating that these systems provide an important source of ionizing photons.

(vi) Due to efficient metal enrichment, Pop III stars are confined to form in the smallest minihaloes at $z > 16$, and their formation rate is larger in the LG than along the MW merger tree. This suggests that traces of Pop III star formation are not confined to the MW and its

<footnote n="18">¹⁸ By definition, the MW merger tree contains only the final MW host halo at $z = 0$.</footnote>
satellites but may be found in external galaxies of the LG, although their detection may be challenging even for the next generation of telescopes.

(vii) We find that a large number of minihaloes having old stellar populations are dragged into the MW or can survive in the local Universe. However, due to the effect of radiative feedback, minihaloes collapsing at $z < z_{reio}$ remain instead dark because they never experienced star formation.

(viii) The low-redshift evolution of all haloes, when followed in stellar mass and metal mass, shows the importance of dynamical effects acting on to progenitors that are being accreted on the major branch of the MW. Events of late mergers, tidal stripping and halo disruptions are found to be relevant in redistributing baryonic properties among halo families, and also prove the capability of GAMESH in tracking the statistical relevance of dynamical effects without accounting for a detailed gas dynamics treatment.

The work that we have presented represents a promising starting point for a more detailed analysis based on sophisticated simulations having a proper treatment of the RT and inhomogeneous metal enrichment with GAMESH.

## ACKNOWLEDGEMENTS

The authors would like to thank Leslie Hunt, Giulia Rodighiero and the anonymous referee for their very constructive comments. LG thanks Michele Ginolfi for his support with observational data. SS was supported by the European Commission through an individual Marie-Skodolowska-Curie Fellowship, project: PRIMORDIAL, 700907. The cosmological simulation and their test simulations for this paper were ran on the UCL facility Grace, and the DiRAC facilities (through the COSMOS and MSSL-Astro consortium) jointly funded by the UK's Science and Technology Facilities Council and the Large Facilities Capital Fund of BIS. We also acknowledge PRACE<fnref n="19" /> for awarding us access to the CEA HPC facility 'CURIE@GENCI'<fnref n="20" /> with the Type B project: High Performance release of the GAMESH pipeline.

The research leading to these results has received funding from the European Research Council under the European Union's Seventh Framework Programme (FP/2007-2013)/ERC Grant Agreement no. 306476.

## REFERENCES

* Abazajian K. et al., 2003, AJ, 126, 2081

* Alam S. et al., 2015, ApJS, 219, 12

* Angulo R. E., Springel V., White S. D. M., Jenkins A., Baugh C. M., Frenk C. S., 2012, MNRAS, 426, 2046

* Behroozi P. S., Wechsler R. H., Conroy C., 2013, ApJ, 770, 57

* Bland-Hawthorn J., Gerhard O., 2016, ARA&A, 54, 529

* Bland-Hawthorn J., Sutherland R., Webster D., 2015, ApJ, 807, 154

* Blanton M. R. et al., 2017, preprint (arXiv:1703.00052)

* Bouwens R. J. et al., 2016, ApJ, 830, 67

* Bovill M. S., Ricotti M., 2009, ApJ, 693, 1859

* Bovy J., Rix H.-W., 2013, ApJ, 779, 115

* Boylan-Kolchin M., Springel V., White S. D. M., Jenkins A., Lemson G., 2009, MNRAS, 398, 1150

* Bromm V., Yoshida N., 2011, ARA&A, 49, 373

* Brook C. B., Kawata D., Scannapieco E., Martel H., Gibson B. K., 2007, ApJ, 661, 10

* Brown T. M. et al., 2014, ApJ, 796, 91

* Caffau E., Ludwig H.-G., Steffen M., Freytag B., Bonifacio P., 2011, Sol. Phys., 268, 255

* Carlesi E. et al., 2016, MNRAS, 458, 900

* Chomiuk L., Povich M. S., 2011, AJ, 142, 197

* Ciardi B., Ferrara A., 2005, Space Sci. Rev., 116, 625

* Ciardi B., Ferrara A., Marri S., Raimondo G., 2001, MNRAS, 324, 381

* Cole S., Aragon-Salamanca A., Frenk C. S., Navarro J. F., Zepf S. E., 1994, MNRAS, 271, 781

* Cole S., Lacey C. G., Baugh C. M., Frenk C. S., 2000, MNRAS, 319, 168

* Colin P., Avila-Reese V., Roca-Fabrega S., Valenzuela O., 2016, ApJ, 829, 98

* Cousin M., Buat V., Boissier S., Bethermin M., Roehlly Y., Génois M., 2016, A&A, 589, A109

* Creasey P., Scannapieco C., Nuza S. E., Yepes G., Gottlöber S., Steinmetz M., 2015, ApJ, 800, L4

* Dayal P., Ferrara A., Dunlop J. S., 2013, MNRAS, 430, 2891

* de Bennassuti M., Schneider R., Valiante R., Salvadori S., 2014, MNRAS, 445, 3039 (DB14)

* de Bennassuti M., Salvadori S., Schneider R., Valiante R., Omukai K., 2017, MNRAS, 465, 926 (DB17)

* De Lucia G., Tornatore L., Frenk C. S., Helmi A., Navarro J. F., White S. D. M., 2014, MNRAS, 445, 970

* Diehl R. et al., 2006, Nature, 439, 45

* Diemand J., Kuhlen M., Madau P., Zemp M., Moore B., Potter D., Stadel J., 2008, Nature, 454, 735

* Fattahi A. et al., 2016, MNRAS, 457, 844

* Flynn C., Holmberg J., Portinari L., Fuchs B., Jahreiß H., 2006, MNRAS, 372, 1149

* Frebel A., Bromm V., 2012, ApJ, 759, 115

* Garrison-Kimmel S., Boylan-Kolchin M., Bullock J. S., Lee K., 2014, MNRAS, 438, 2578

* Geen S., Slyz A., Devriendt J., 2013, MNRAS, 429, 633

* Graziani L., Maselli A., Ciardi B., 2013, MNRAS, 431, 722

* Graziani L., Salvadori S., Schneider R., Kawata D., de Bennassuti M., Maselli A., 2015, MNRAS, 449, 3137 (LG15)

* Griffen B. F., Dooley G. A., Ji A. P., O'Shea B. W., Gómez F. A., Frebel A., 2016a, preprint (arXiv:1611.00759)

* Griffen B. F., Ji A. P., Dooley G. A., Gómez F. A., Vogelsberger M., O'Shea B. W., Frebel A., 2016b, ApJ, 818, 10

* Grogin N. A. et al., 2011, ApJS, 197, 35

* Guo Q., White S., Li C., Boylan-Kolchin M., 2010, MNRAS, 404, 1111

* Guo Y. et al., 2016, ApJ, 822, 103

* Hahn O., Abel T., 2011, MNRAS, 415, 2101

* Hartwig T., Bromm V., Klessen R. S., Glover S. C. O., 2015, MNRAS, 447, 3892

* Henriques B. M. B., White S. D. M., Thomas P. A., Angulo R., Guo Q., Lemson G., Springel V., Overzier R., 2015, MNRAS, 451, 2663

* Hirschmann M., Naab T., Somerville R. S., Burkert A., Oser L., 2012, MNRAS, 419, 3200

* Hirschmann M. et al., 2013, MNRAS, 436, 2929

* Hopkins P. F., Kereš D., Óñorbe J., Faucher-Giguère C.-A., Quataert E., Murray N., Bullock J. S., 2014, MNRAS, 445, 581

* Huchra J., Davis M., Latham D., Tonry J., 1983, ApJS, 52, 89

* Hunt L. et al., 2012, MNRAS, 427, 906

* Hunt L., Dayal P., Magrini L., Ferrara A., 2016a, MNRAS, 463, 2002

* Hunt L., Dayal P., Magrini L., Ferrara A., 2016b, MNRAS, 463, 2020

* Ibata R., Martin N. F., Irwin M., Chapman S., Ferguson A. M. N., Lewis G. F., McConnachie A. W., 2007, ApJ, 671, 1591

* Kawata D., Gibson B. K., 2003, MNRAS, 340, 908

* Kawata D., Okamoto T., Gibson B. K., Barnes D. J., Cen R., 2013, MNRAS, 428, 1968

* Kennicutt R. C., Jr, 1998, ARA&A, 36, 189

* Kennicutt R. C., Evans N. J., 2012, ARA&A, 50, 531

* Kewley L. J., Ellison S. L., 2008, ApJ, 681, 1183

* Kim J.-h. et al., 2014, ApJS, 210, 14

* Kirby E. N., Simon J. D., Geha M., Guhathakurta P., Frebel A., 2008, ApJ, 685, L43

<footnote n="19">¹⁹ http://www.prace-ri.eu/</footnote>

<footnote n="20">²⁰ http://www-hpc.cea.fr/en/complexe/tgcc-curie.htm</footnote>
Klypin A. A., Trujillo-Gomez S., Primack J., 2011, ApJ, 740, 102
Knebe A. et al., 2015, MNRAS, 451, 4029
Komatsu E. et al., 2009, ApJS, 180, 330
Komiya Y., Suda T., Minaguchi H., Shigeyama T., Aoki W., Fujimoto
M. Y., 2007, ApJ, 658, 367
Kubryk M., Prantzos N., Athanassoula E., 2015, A&A, 580, A126
Larson R. B., 1998, MNRAS, 301, 569
Larson D. et al., 2011, ApJS, 192, 16
Licquia T. C., Newman J. A., 2015, ApJ, 806, 96
McConnachie A. W., 2012, AJ, 144, 4
Mckee C. F., Ostriker E. C., 2007, ARA&A, 45, 565
McKee C. F., Williams J. P., 1997, ApJ, 476, 144
McMillan P. J., 2011, MNRAS, 414, 2446
Madden S. C. et al., 2013, PASP, 125, 600
Madden S. C. et al., 2014, PASP, 126, 1079
Maiolino R. et al., 2008, A&A, 488, 488
Mannucci F. et al., 2009, MNRAS, 398, 1915
Mannucci F., Cresci G., Maiolino R., Marconi A., Gnerucci A., 2010,
MNRAS, 408, 2115
Maselli A., Ferrara A., Ciardi B., 2003, MNRAS, 345, 379
Maselli A., Ciardi B., Kanekar A., 2009, MNRAS, 393, 171
Ménard B., Scranton R., Fukugita M., Richards G., 2010, MNRAS, 405,
1025
Mo H., van den Bosch F. C., White S., 2010, Galaxy Formation and Evolu-
tion. Cambridge Univ. Press, Cambridge
Monelli M. et al., 2016, ApJ, 819, 147
Muñoz J. A., Madau P., Loeb A., Diemand J., 2009, MNRAS, 400, 1593
Navarro J. F., Frenk C. S., White S. D. M., 1996, ApJ, 462, 563
Navarro J. F. et al., 2010, MNRAS, 402, 21
Nuza S. E., Parisi F., Scannapieco C., Richter P., Gottlöber S., Steinmetz
M., 2014, MNRAS, 441, 2593
Ocvirk P. et al., 2014, ApJ, 794, 20
Papovich C. et al., 2015, ApJ, 803, 26
Patel S. G. et al., 2013, ApJ, 778, 115
Peeples M. S., Werk J. K., Tumlinson J., Oppenheimer B. D., Prochask
J. X., Katz N., Weinberg D. H., 2014, ApJ, 786, 54
Pezzulli E., Valiante R., Schneider R., 2016, MNRAS, 458, 3047
Planck Collaboration XVI, 2014, A&A, 571, A16
Powell L. C., Slyz A., Devriendt J., 2011, MNRAS, 414, 3671
Press W. H., Schechter P., 1974, ApJ, 187, 425
Robitaille T. P., Whitney B. A., 2010, ApJ, 710, L11
Salvadori S., Ferrara A., 2009, MNRAS, 395, L6
Salvadori S., Ferrara A., 2012, MNRAS, 421, L29
Salvadori S., Schneider R., Ferrara A., 2007, MNRAS, 381, 647
Salvadori S., Ferrara A., Schneider R., Scannapieco E., Kawata D., 2010a,
MNRAS, 401, L5
Salvadori S., Dayal P., Ferrara A., 2010b, MNRAS, 407, L1
Salvadori S., Tolstoy E., Ferrara A., Zaroubi S., 2014, MNRAS, 437,
L26
Salvadori S., Skúladóttir Á., Tolstoy E., 2015, MNRAS, 454, 1320
Sarmento R., Scannapieco E., Pan L., 2017, ApJ, 834, 23
Sawala T. et al., 2016, MNRAS, 457, 1931
Scannapieco E., Schneider R., Ferrara A., 2003, ApJ, 589, 35
Scannapieco E., Kawata D., Brook C. B., Schneider R., Ferrara A., Gibson
B. K., 2006, ApJ, 653, 285
Scannapieco C., White S. D. M., Springel V., Tissera P. B., 2011, MNRAS,
417, 154
Scannapieco C. et al., 2012, MNRAS, 423, 1726
Scannapieco C., Creasey P., Nuza S. E., Yepes G., Gottlöber S., Steinmetz
M., 2015, A&A, 577, A3
Schaye J. et al., 2015, MNRAS, 446, 521
Schreiber C. et al., 2015, A&A, 575, A74
Simon J. D., Geha M., 2007, ApJ, 670, 313
Skibba R. A. et al., 2012, ApJ, 761, 42
Smith L. F., Biermann P., Mezger P. G., 1978, A&A, 66, 65
Springel V. et al., 2005, Nature, 435, 629
Springel V., Frenk C. S., White S. D. M., 2006, Nature, 440, 1137
Springel V. et al., 2008, MNRAS, 391, 1685

Terrazas B. A., Bell E. F., Henriques B. M. B., White S. D. M., 2016, MNRAS, 459, 1929
Tolstoy E., Hill V., Tosi M., 2009, ARA&A, 47, 371
Tumlinson J., 2006, ApJ, 641, 1
Tumlinson J., 2010, ApJ, 708, 1398
Valiante R., Schneider R., Salvadori S., Bianchi S., 2011, MNRAS, 416, 1916
Valiante R., Schneider R., Salvadori S., Gallerani S., 2014, MNRAS, 444, 2442
Valiante R., Schneider R., Volonteri M., Omukai K., 2016, MNRAS, 457, 3356
van Dokkum P. G. et al., 2013, ApJ, 771, L35
Vincenzo F., Matteucci F., Vattakunnel S., Lanfranchi G. A., 2014, MNRAS, 441, 2815
Vogelsberger M. et al., 2014, Nature, 509, 177
Wang J. et al., 2011, MNRAS, 413, 1373
Wang W., Han J., Cooper A. P., Cole S., Frenk C., Lowing B., 2015, MNRAS, 453, 377
Wetzel A. R., Deason A. J., Garrison-Kimmel S., 2015, ApJ, 807, 49
White S. D. M., Frenk C. S., 1991, ApJ, 379, 52
 Yepes G., Gottlöber S., Hoffman Y., 2014, New Astron. Rev., 58, 1

APPENDIX A: MW HALO PROPERTIES
AND COMPARISON WITH OTHER ZOOM-IN
SIMULATIONS

In this appendix we summarize the structural properties of the MW and compare with similar haloes found in the *AQUARIUS*, *ELVIS*, *CATERPILLAR* and *APOSTLE* simulations. Tables A1 and A2 collect haloes having $M_{vir} \sim 1.5-1.9 \times 10^{12} M_{\odot}$ and summarize their structural properties. The Hamilton/Hall haloes of the *ELVIS* catalogue and Cat-8/Cat-12 in the *CATERPILLAR* have been selected as comparison targets from the available statistical samples, while the properties of Aq-A and Aq-C will be discussed later, when comparing their accretion history to the MW one. We also included halo Ap-8 from the *APOSTLE* project even if this is a hydrodynamical simulation and many structural properties of Ap-8 are not documented. Its satellite distribution, on the other hand, has been extensively studied in the literature and will allow future dedicated comparisons (Mancini et al., in preparation).

<table><caption>**Table A1.** Summary of reference zoom-in simulations adopted for comparison with the MW halo found in *GAMESH: AQUARIUS* (Springel et al. 2008), *ELVIS* (Garrison-Kimmel et al. 2014), *CATERPILLAR* (Griffen et al. 2016b,a), *APOSTLE* (Fattahi et al. 2016; Sawala et al. 2016). In this table 'Halo ID' is the unique name of the halo as found in the cited literature, while the column named 'Pair' indicates if the halo is found in a binary configuration with a M31 analogue or it is an isolated one. The *Planck* cosmology refers to Planck Collaboration XVI (2014) while *WMAP-5* to Komatsu et al. (2009) and *WMAP-7* to Larson et al. (2011).</caption><thead><tr><th>Simulation</th><th>Halo ID</th><th>Pair</th><th>Cosmology</th></tr></thead><tbody><tr><td>GAMESH</td><td>MW</td><td>No</td><td>Planck</td></tr><tr><td>AQUARIUS</td><td>Aq-A</td><td>No</td><td>WMAP-5</td></tr><tr><td>AQUARIUS</td><td>Aq-C</td><td>No</td><td>WMAP-5</td></tr><tr><td>ELVIS</td><td>Hamilton</td><td>Yes</td><td>WMAP-7</td></tr><tr><td>ELVIS</td><td>iHamilton</td><td>No</td><td>WMAP-7</td></tr><tr><td>ELVIS</td><td>Hall</td><td>Yes</td><td>WMAP-7</td></tr><tr><td>ELVIS</td><td>iHall</td><td>No</td><td>WMAP-7</td></tr><tr><td>ELVIS</td><td>iHall HiRes</td><td>No</td><td>WMAP-7</td></tr><tr><td>CATERPILLAR</td><td>Cat-8/LX13</td><td>No</td><td>Planck</td></tr><tr><td>CATERPILLAR</td><td>Cat-8/LX14</td><td>No</td><td>Planck</td></tr><tr><td>CATERPILLAR</td><td>Cat-12/LX13</td><td>No</td><td>Planck</td></tr><tr><td>CATERPILLAR</td><td>Cat-12/LX14</td><td>No</td><td>Planck</td></tr><tr><td>APOSTLE/DOVE</td><td>Ap-8</td><td>Yes</td><td>WMAP-7</td></tr></tbody></table>
<table><caption>**Table A2.** Summary of MW-like halo properties taken for comparison with *GAMESH* (see Table A1). Here 'Halo ID' is the unique name of the halo, $m_p$ is the DM particle mass adopted for the highest resolution run, $M_{vir}$ and $R_{vir}$ are virial mass and radius respectively, $R_s$ is the maximum radius of the sphere non-contaminated by lower resolution particles. (Note that in *GAMESH* this value refers to half of the side length of the cube contaminated only by high-resolution particles.) $z_{0.5}$ is the halo formation redshift (see text for more details). The concentration parameter $c$ is calculated following Navarro et al. (1996, 2010), while $V_{max}$ and $R_{max}$ are computed from the rotation curve of the halo (see Garrison-Kimmel et al. 2014 for a definition). Finally, note that some values are not found in the reference literature of the *AQUARIUS*, *ELVIS* and *APOSTLE* simulations.</caption><thead><tr><th>Halo ID</th><th><i>m</i><sub>p</sub> (M<sub>☉</sub>)</th><th><i>M</i><sub>vir</sub> (M<sub>☉</sub>)</th><th><i>R</i><sub>vir</sub> (kpc)</th><th><i>R</i><sub>s</sub> (Mpc)</th><th><i>z</i><sub>0.5</sub></th><th><i>c</i></th><th><i>V</i><sub>max</sub> (km s<sup>-1</sup>)</th><th><i>R</i><sub>max</sub> (kpc)</th></tr></thead><tbody><tr><td>MW</td><td>3.38 × 10<sup>5</sup></td><td>1.72 × 10<sup>12</sup></td><td>317</td><td>2.0</td><td>1.46</td><td>12.7</td><td>198.19</td><td>56.5</td></tr><tr><td>Aq-A-4</td><td>3.93 × 10<sup>5</sup></td><td>1.84 × 10<sup>12</sup></td><td>246</td><td>–</td><td>1.93</td><td>16.21</td><td>209.24</td><td>28.2</td></tr><tr><td>Aq-A-2</td><td>1.37 × 10<sup>4</sup></td><td>1.84 × 10<sup>12</sup></td><td>246</td><td>–</td><td>1.93</td><td>16.19</td><td>208.49</td><td>28.1</td></tr><tr><td>Aq-C-4</td><td>3.21 × 10<sup>5</sup></td><td>1.79 × 10<sup>12</sup></td><td>244</td><td>–</td><td>2.23</td><td>14.84</td><td>223.20</td><td>33.6</td></tr><tr><td>Aq-C-2</td><td>1.40 × 10<sup>4</sup></td><td>1.77 × 10<sup>12</sup></td><td>243</td><td>–</td><td>2.23</td><td>15.21</td><td>222.40</td><td>32.5</td></tr><tr><td>Hamilton</td><td>1.90 × 10<sup>5</sup></td><td>1.76 × 10<sup>12</sup></td><td>315</td><td>1.39</td><td>1.47</td><td>9.9</td><td>197</td><td>–</td></tr><tr><td>iHamilton</td><td>1.90 × 10<sup>5</sup></td><td>1.86 × 10<sup>12</sup></td><td>321</td><td>1.55</td><td>2.11</td><td>14.2</td><td>203</td><td>–</td></tr><tr><td>Hall</td><td>1.90 × 10<sup>5</sup></td><td>1.52 × 10<sup>12</sup></td><td>299</td><td>1.35</td><td>1.04</td><td>10.3</td><td>180</td><td>–</td></tr><tr><td>iHall</td><td>1.90 × 10<sup>5</sup></td><td>1.71 × 10<sup>12</sup></td><td>300</td><td>1.59</td><td>1.13</td><td>6.0</td><td>172</td><td>–</td></tr><tr><td>iHall HiRes</td><td>2.35 × 10<sup>4</sup></td><td>1.67 × 10<sup>12</sup></td><td>309</td><td>1.59</td><td>1.07</td><td>5.8</td><td>167</td><td>–</td></tr><tr><td>Cat-8/LX13</td><td>2.39 × 10<sup>5</sup></td><td>1.70 × 10<sup>12</sup></td><td>315</td><td>1.55</td><td>1.52</td><td>13.3</td><td>197.64</td><td>39.81</td></tr><tr><td>Cat-8/LX14</td><td>2.98 × 10<sup>4</sup></td><td>1.70 × 10<sup>12</sup></td><td>315</td><td>1.54</td><td>1.52</td><td>13.5</td><td>198.56</td><td>40.82</td></tr><tr><td>Cat-12/LX13</td><td>2.39 × 10<sup>5</sup></td><td>1.77 × 10<sup>12</sup></td><td>319</td><td>1.239</td><td>1.37</td><td>11.7</td><td>191.32</td><td>49.44</td></tr><tr><td>Cat-12/LX14</td><td>2.98 × 10<sup>4</sup></td><td>1.76 × 10<sup>12</sup></td><td>319</td><td>1.162</td><td>1.37</td><td>11.4</td><td>191.30</td><td>52.72</td></tr><tr><td>Ap-8</td><td>5.0 × 10<sup>5</sup>/8.8 × 10<sup>6</sup></td><td>1.72 × 10<sup>12</sup></td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td></tr></tbody></table>

From Table A2 we infer that our simulation adopts a mass resolution and cosmology compatible with level 13 (Haloes Cat-8/LX13 and Cat-12/LX13) in the CATERPILLAR project, while the ELVIS simulation provides haloes better resolved by a factor of 1.8. The latter two simulations have released halo catalogues with a mass resolution increased roughly by one order of magnitude: iHall HighRes and Cat-8/LX14, Cat-12/LX14. As for the adopted cosmology, our simulation and the CATERPILLAR suite adopt the same Planck parameters while ELVIS and APOSTLE rely on 7-year Wilkinson Microwave Anisotropy Probe (WMAP-7) measurements. Finally, ELVIS has simulated both isolated (iHamilton/iHall) and paired haloes (with an M31 analogue at a distance of ~800 kpc), while the MW halo and the haloes identified in the CATERPILLAR sample are isolated.

By checking the values in Table A2 it is evident that the MW halo structural properties, such as its virial radius $R_{vir}$, concentration parameter $c$ (see Navarro et al. 1996, 2010 for a definition) and maximum circular velocity $V_{max}$ (Garrison-Kimmel et al. 2014) are compatible with the scatter in the Cat-8/Cat-12 at level 13 and iHamilton/iHall. Interestingly, we note that the larger scatter in halo structural properties($c$, $R_{vir}$, $R_{max}$) can be linked to their peculiar accretion histories (see Section 2.4 for more details).

The ability of these simulations to resolve the halo environment is quantified by $R_s$, the size of the 'non-contaminated region'. In our simulation, a cubic volume of 4 cMpc side length is uncontaminated, which is larger than spherical regions with $R_s = 1.1-1.5$ cMpc surrounding the other haloes. This rich (but still computationally affordable) dynamical information will allow future studies to adopt accurate RT feedback and/or particle-tagging techniques to model inhomogeneous enrichment and to trace the distribution of stellar populations inside the MW halo.

## APPENDIX B: COMPARISON WITH OTHER SAM PREDICTIONS

Here we briefly compare our predictions on the redshift evolution of SFR, $M_*$ and $M_Z$, with independent SAMs coupled to DM simulations. As often discussed in many SAM comparison projects (Hirschmann et al. 2012; Knebe et al. 2015) a clear comparison

should be based on a set of similar or, at least, controlled physical assumptions. This is largely beyond the scope of this paper because our model is based on single DM halo simulation and does not provide any statistical prediction across the natural diversity in the growth of DM haloes (Navarro et al. 2010; Terrazas et al. 2016). Despite this intrinsic limitation, when compared with predictions based on MW-like haloes having similar mass, structural properties and smooth accretion history, *GAMESH* results are in broad agreement with the mean SFR and $M_*$ shown in Hirschmann et al. (2013) (see top left-hand panels of figs 2 and 5 relative to $12 < \log M < 13$) or with the FULL case of fig. 6 in Hirschmann et al. (2012). Note that these models are not tuned to reproduce the MW properties at $z=0$ (as imposed in *GAMESH*) and then their final values can be sensitively different.

Once *GAMESH* is compared with predictions based on Aq-A-5, Aq-C-5 (Scannapieco et al. 2011), our model agrees with the $M_*$ evolution even if the intrinsic differences in the SAM assumptions and parameter tuning cannot simultaneously guarantee a similar trend in SFR. This is also the case for Aq-A-3 discussed in De Lucia et al. (2014).

An excellent agreement is found instead with the model $m_3$ of Cousin et al. (2016), mainly because it adopts a similar description of star formation and metal enrichment. These authors use in fact an identical prescription to compute the SFR and adopt an efficiency $\epsilon_*$ = 0.02 on a set of 56 MW-like haloes having $3 < M_*/10^{10} M_{\odot} < 7$ and $0.7 < M_{DM}/10^{12} M_{\odot} < 3$ (in their jargon MW-sister and MW-cousins haloes). They also find SFR($z=0$) = (1-5) $M_{\odot} yr^{-1}$. Despite their disc modelling is not implemented in *GAMESH*, our MW halo can be certainly classified as a MW-cousin and the comparison of model predictions for the SFR is simply straightforward (compare fig. 6 with fig. 9 in their paper).

This paper has been typeset from a TeX/LATeX file prepared by the author.