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

TL;DR

This study models the evolution of Milky Way-like galaxies from early cosmic times to the present, predicting properties of their progenitors and analyzing the roles of different halo types in galaxy formation and reionization.

Contribution

It applies a novel simulation approach to trace MW progenitors from high redshift, revealing the significance of mini-halos and dynamical effects in galaxy evolution and reionization.

Findings

MW progenitors follow observed relations up to z=4

Mini-halos contribute significantly to early star formation

Dynamical effects redistribute baryonic matter at low redshift

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 (Graziani et al. 2015) 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 halos, at z > 6 the contribution of star forming mini-halos 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 mini-halos 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.

Tables

Table 1. Table 1: 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 ; Griffen et al., 2016a ) , 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 et al. ( 2014 ) while WMAP-5 to Komatsu et al. ( 2009 ) and WMAP-7 to Larson et al. ( 2011 ) .

Simulation	Halo ID	Pair	Cosmology
GAMESH	MW	no	Planck
AQUARIUS	Aq-A	no	WMAP-5
AQUARIUS	Aq-C	no	WMAP-5
ELVIS	Hamilton	yes	WMAP-7
ELVIS	iHamilton	no	WMAP-7
ELVIS	Hall	yes	WMAP-7
ELVIS	iHall	no	WMAP-7
ELVIS	iHall HiRes	no	WMAP-7
CATERPILLAR	Cat-8/LX13	no	Planck
CATERPILLAR	Cat-8/LX14	no	Planck
CATERPILLAR	Cat-12/LX13	no	Planck
CATERPILLAR	Cat-12/LX14	no	Planck
APOSTLE/DOVE	Ap-8	yes	WMAP-7

Table 2. Table 2: Summary of Milky Way-like halo properties taken for comparison with GAMESH (see Table  1 ). 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 vir {}_{\texttt{vir}} , R vir vir {}_{\texttt{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 subscript 𝑧 0.5 z_{0.5} is the halo formation redshift (see text for more details). The concentration parameter c 𝑐 c is calculated following Navarro et al. 1996 , 2010 , while V max max {}_{\texttt{max}} and R max max {}_{\texttt{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.

Halo ID	mp $\left[M_{\odot }\right]$	M${}_{\text{vir}}$ $\left[M_{\odot }\right]$	R${}_{\text{vir}}$ [kpc]	Rs [Mpc]	$z_{0.5}$	$c$	V${}_{\text{max}}$ [km/s]	R${}_{\text{max}}$ [kpc]
MW	$3.38\times {10}^5$	$1.72\times {10}^{12}$	317	2.0	1.46	12.7	198.19	56.5
Aq-A-4	3.93$\times {10}^5$	$1.84\times {10}^{12}$	246	-	1.93	16.21	209.24	28.2
Aq-A-2	1.37$\times {10}^4$	$1.84\times {10}^{12}$	246	-	1.93	16.19	208.49	28.1
Aq-C-4	3.21$\times {10}^5$	$1.79\times {10}^{12}$	244	-	2.23	14.84	223.20	33.6
Aq-C-2	1.40$\times {10}^4$	$1.77\times {10}^{12}$	243	-	2.23	15.21	222.40	32.5
Hamilton	1.90$\times {10}^5$	$1.76\times {10}^{12}$	315	1.39	1.47	9.9	197	-
iHamilton	1.90$\times {10}^5$	$1.86\times {10}^{12}$	321	1.55	2.11	14.2	203	-
Hall	1.90$\times {10}^5$	$1.52\times {10}^{12}$	299	1.35	1.04	10.3	180	-
iHall	1.90$\times {10}^5$	$1.71\times {10}^{12}$	300	1.59	1.13	6.0	172	-
iHall HiRes	2.35$\times {10}^4$	$1.67\times {10}^{12}$	309	1.59	1.07	5.8	167	-
Cat-8/LX13	2.39$\times {10}^5$	$1.70\times {10}^{12}$	315	1.55	1.52	13.3	197.64	39.81
Cat-8/LX14	2.98$\times {10}^4$	$1.70\times {10}^{12}$	315	1.54	1.52	13.5	198.56	40.82
Cat-12/LX13	2.39$\times {10}^5$	$1.77\times {10}^{12}$	319	1.239	1.37	11.7	191.32	49.44
Cat-12/LX14	2.98$\times {10}^4$	$1.76\times {10}^{12}$	319	1.162	1.37	11.4	191.30	52.72
Ap-8	5.0$\times {10}^5$/8.8$\times {10}^6$	$1.72\times {10}^{12}$	-	-	-	-	-	-

Full text

††pagerange: The history of the dark and luminous side of Milky Way-like progenitors–B††pubyear: 2016

1 Introduction

Modern cosmological models (White & Frenk, 1991; Cole et al., 1994; Cole et al., 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 et al., 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 at the center of any model of galaxy formation and evolution.

The possibility of resolving stars, both in the Milky Way 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; Brook et al., 2007; Komiya et al., 2007; Tumlinson, 2010; Salvadori et al., 2010a; Frebel & Bromm, 2012; de Bennassuti et al., 2014; Hartwig et al., 2015; Salvadori et al., 2015; de Bennassuti et al., 2016).

The great number of kinematic and chemical tracers available in the Milky Way (see Bland-Hawthorn & Gerhard 2016 for a recent review) have 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 et al. 2009; McConnachie 2012; Madden et al. 2013; Monelli et al. 2016). Ultra-faint 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 (Salvadori & Ferrara, 2009; Bovill & Ricotti, 2009; Muñoz et al., 2009; Vincenzo et al., 2014; Bland-Hawthorn et al., 2015; Salvadori et al., 2015). Given this wealth of data, theoretical models have attempted to understand how the star formation history of the Milky Way and its dwarf companions (Kennicutt & Evans, 2012) is affected by large scale processes like cosmic reionization and metal enrichment (Ocvirk et al., 2014; Graziani et al., 2015; Wetzel et al., 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 Millenium suite111http://wwwmpa.mpa-garching.mpg.de/millennium/ (Springel et al., 2005, 2006; Boylan-Kolchin et al., 2009; Angulo et al., 2012) or the MultiDark and Bolshoi runs222https://www.cosmosim.org/ (Klypin et al., 2011) have traditionally provided the theoretical framework on top of which semi-analytic models interpret global observed quantities (see for instance Knebe et al. 2015). Recently, two hydrodynamical projects (Schaye et al., 2015; Vogelsberger et al., 2014) including detailed feedback on sub-grid scales have been able to reach an unprecedented realism in reproducing morphological and structural galaxy properties.

In the past years many dark matter simulations of MW-like halos 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 (Yepes et al., 2014; Nuza et al., 2014) and the ELVIS simulation suite (Garrison-Kimmel et al., 2014). They have been extensively used to study the structural properties of dark matter halos, the statistics of their satellites, as well as to correctly constrain the initial conditions leading to the dynamical configuration of our Local Universe, having a MW-M31 galaxy pair (Carlesi et al., 2016). Halo assembly histories and the role of mergers 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, semi-analytic models (SAM) 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; Griffen et al. 2016b).

Despite a long series of investigations on how to implement mechanical and thermal feedback in different hydro-dynamical schemes (Powell et al., 2011; Geen et al., 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 favor 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 photo-ionization modeling is preferable (Graziani et al., 2013). Finally, some observable quantities can be simply more sensitive to photo-ionization 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 ("GAlaxy MErger Tree and Evolution", Salvadori et al. 2007, 2010a) has been considerably extended (see de Bennassuti et al. 2014, hereafter DB14) and self-consistently coupled with the radiative transfer code CRASH ("Cosmological RAdiative transfer Scheme for Hydrodynamics", Ciardi et al. 2001; Maselli et al. 2003, 2009; Graziani et al. 2013) creating GAMESH, the first implementation of a full radiative feedback treatment in a semi-analitic model on top of a DM only simulation. GAMESH is then capable to span the entire formation process of a Milky Way-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 Milky Way halo and its satellites, the formation and coalescence environments of massive black hole binaries, the molecular and dust content of Milky-Way 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, Milky Way-type halo in a cosmic cube of about 4 cMpc side length. The application of the semi-analytic model 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 Milky Way 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 Milky Way progenitor galaxies and we critically compare them with observations of Milky-Way 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; Hunt et al., 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 Milky Way and its Local Group. Furthermore, it allows 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 catalog, 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), APOSTLE (Fattahi et al., 2016; Sawala et al., 2016) (see Table 1 in Appendix A for more details).

Section 3 describes the tuning of GAMESH on a set of observed Milky Way properties. Section 4 focuses on analyzing the properties of MW progenitors, while Section 5 analyses their evolution as hosted in mini-halos or Lyman $\alpha$-cooling halos. 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 a 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 as well as the statistics of various halo populations in a surrounding volume of 4 cMpc side length. A careful comparison of the MW properties with similar halos taken from independent simulations both in single and paired configurations can be found in Appendix A, where the halo properties are summarized in Table 2 for an easier comparison.

2.1 GCD+ and initial conditions

The N-body cosmological simulation of a MW-sized halo has been performed with GCD+ (Kawata & Gibson, 2003; Kawata et al., 2013) with a $\beta$-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 (Planck Collaboration et al., 2014) ( $\Omega_{0}=0.32,\Lambda_{0}=0.78,\Omega_{b}=0.049$ and $h=0.67$) to simulate a volume of $83.5^{3}$ Mpc3. In this volume we identified a Milky Way-sized halo and we created the initial conditions for a a zoom-in simulation. The final run consists of a total of 62421192 particles, 55012200 of which in the highest resolution region having particles with mass of $3.4\times 10^{5}$ M*⊙. The virial mass of the resulting MW halo is $1.7\times 10^{12}$ M⊙*.

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 radiative transfer is adopted.

The final output of the simulation provides a list of collapsed halo objects as well as the projection onto grids of 5123 cells/side of the DM distribution found in the 4 cMpc cosmic volume centered 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 halos 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_{\rm vir}$, mass $M_{\rm vir}$ and radius $R_{\rm vir}$), as well the dynamical variables (position and velocity) of its center of mass with respect to the central MW halo. Besides the list of DM halos present in each snapshots, we also stored position and velocities of all particles resolving them. This information allows to study the internal structure of the most massive halos found in the LG (DM profiles, angular momentum, internal motion, over-density structures, etc..), once they have a sufficient number of particles to reliably compute these quantities. A detailed dynamical study of the MW satellite and sub-satellite properties is still in preparation (see Mancini et al., in prep.). 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_{\texttt{DM}}\sim 3\times 10^{12}M_{\odot}$ distributed in $2458$ halos. Among these, 2 halos have DM mass $M\sim 10^{11}$ M*⊙* (M32, M33 or LMC-type halos, see Table 1 in Guo et al. 2010 and references therein), 14 have $10^{10}\lesssim M<10^{11}$ M*⊙, and 98 have $10^{9}\lesssim M<10^{10}$ M⊙* (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_{\rm M31}\gtrsim M_{\rm MW}$ by recent estimates (see Ibata et al. 2007 and references therein.). Other M31-sized halos 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 halos 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 OpenMP333www.openmp.org 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}]$444Note that once a particle belongs to many halos, its multiplicity is also computed and stored to exactly account for the particle contribution in mass transfer processes across the merger tree steps.. 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 re-centering adopted for each snapshot to define the LG volume in the simulation data, few pIDs are sporadically not found in the LG volume because 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 halos 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 halos with their environments and to exactly follow all their dynamical processes regulating the accretion of dark matter halos: mergers, tidal stripping and halo disruptions. All these events can then be classified and analyzed 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 IGM, by scaling with their relative dark mass contribution. While this is the approach used in the 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 halos with 100 particles, our IGM could contain additional population of halos with unresolved masses.

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

• •

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

• •

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

• •

Halo stripping: when a halo loses part of its mass by tidal interactions with nearby halos.

• •

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, Milky Way-like, dark matter halo (hereafter MW) found at the center 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 Local Group.

Figure 1 shows the total mass of the MW halo 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 halos providing particles collapsing onto the MW by the successive snapshot. As 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_{\texttt{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 halos. 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_{\texttt{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_{\texttt{MM}}(z)=M_{\texttt{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 dark matter simulations described in Behroozi et al. (2013, see in particular their Figure 19), and relative to halos of similar mass. A further comparison with Table 2 shows that the MW halo has a formation redshift z0.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 different accretion and dynamical history. Also note that the same time obtained from the blue line (i.e. accounting for all progenitors, see Navarro et al. 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_{\texttt{DM}}(a)/M_{\texttt{MW}}(a=1)$, in the bottom left corner of this figure, we show the same history in these units for a more straightforward comparison with other simulations (see for example Figure 9 in Griffen et al. 2016b).

The resulting MW mass found at $z=0$ ($M_{\texttt{MW}}\sim 1.7\times 10^{12}M_{\odot}$) is finally compared with the scatter in mass of MW-size halos555The mass ranges have been computed by extracting min/max mass values in tables of relative papers. Halos 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. found in recent simulations targeting the Local Group, 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 hydro-dynamical APOSTLE simulation (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 Figure 1 and cross-check with the visual picture provided by Figure 2, which shows the redshift evolution of the LG in a series of slice cuts intercepting the central MW galaxy. Here the DM density map, is shown as gradient from white (collapsed regions) to black (voids)666Note that the equivalent gas number density is obtained by scaling the DM mass by the universal baryon fraction.

It is immediately evident that while at high redshifts the mass of both MW and LG have a similar evolution777This is mainly because all the halos collapse first along a filament at the center of the box., below $t=11$ Gyr ($z\sim 2$) many structures not belonging to the MW merger tree, start collapsing in the entire volume or enter the domain from the larger scale888This is a $8$ cMpc cube assumed to gravitationally constrain the LG domain and contains structures described by lower resolution particles.. The evolution at high redshifts proceeds by assembling halos 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\sim 3$, the central halo dynamically dominates the LG region and continues to drag material entering from larger scales: around $z\sim 2$, an external filament not previously visible within $4$ cMpc provides halos to the central galaxy. This is a clear hint that galaxy formation is a multi-scale 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 attractor999An animation can be found in the article online resource files..

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 Figure 3 where we show the percentage of mass increase relative to the final MW mass ($\Delta M(t)/M_{\texttt{MW}}(z=0)$), as a function of the lookback time $t$ (solid red line). It is immediately evident that across cosmic times, the Milky Way halo grows by means of a smooth and continuous assembly of matter spaced out by many violent accretion events, each of which provides a $\sim$3% contribution to the final mass (see for example 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%, 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 Lyman $\alpha$ (Ly$\alpha$-) cooling halos ($T_{\rm vir}>2\times 10^{4}$ K) found in the MW merger tree101010Note 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$.. 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% to the final mass. Note that minor mergers are continuously found between the MW and small halos or between Lyman $\alpha$ cooling halos and mini-halos ($T_{\rm 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 % to the final dark matter 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 Millenium-II simulation (Boylan-Kolchin et al., 2009), the growth of dark matter halos can be dominated by mergers or by a smooth accretion of diffuse matter. Wang et al. (2011) consider 6 halos of the Aquarius simulation (Springel et al., 2008) which target a Milky Way-like halo mass at $z=0$ ($1-2\times 10^{12}$ M*⊙*) and study their accretion mode and its impact on the internal structure and age distributions of particles in the final Halo. The authors claim that by averaging over the 6 halos, smooth accretion can provide a relevant contribution to the final mass (roughly 30-40%). 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%, in agreement with the estimates discussed above for MW111111Note that in this work we adopted a ratio (1:4) to identify a major merger event.. 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 mini-halos and Lyman $\alpha$-cooling halos. 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, and de Bennassuti et al. 2016, hereafter DB16).

Figure 4 shows the redshift evolution of the number of mini-halos (thin lines) and Ly$\alpha$-cooling halos (thick lines) along the merger tree of the MW (blue solid lines) and in the LG (black dotted lines). First note the plethora of mini-halos predicted by the DM simulation around the MW at $z=0$ (N${}_{\rm Mini}\sim 550$ within $2\times R_{\rm vir}$). In fact, during the last step of its mass assembly, less than 30% of the entire mini-halo population in the LG volume is embedded in the final MW halo; these halos 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 halos, their number is about one order of magnitude lower ($\sim 3$%), with very few objects (N${}_{\rm 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 halos 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 halos and by assuming an instant reionization at $z_{\rm reio}=6$121212Alternative, redshift modulated expressions can be found in Salvadori & Ferrara 2009, 2012, DB14, LG15, DB16.. 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; Salvadori et al., 2014; de Bennassuti et al., 2014, 2016) or adapt the original code to investigate the formation of quasars (Valiante et al., 2011, 2014, 2016a; Pezzulli et al., 2016), the pure semi-analytic scheme of GAMESH is based on the following simplifying assumptions:

• •

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

• •

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

• •

Stellar evolution is followed assuming the Instantaneous Recycling Approximation (IRA). When the gas metallicity is $Z<Z_{\rm cr}=10^{-4}\,Z_{\odot}$, Population III stars are formed with a constant mass of 200 $M_{\odot}$. Above the critical metallicity, Population 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_{\rm ch}=0.35M_{\odot}$.

• •

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.

• •

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

• •

When $z\leq z_{\rm reio}$ star formation can only occur in galaxies with $T_{\rm vir}>2\times 10^{4}$ K, to account for the effects of photo-heating and photo-evaporation (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 DB16) 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 Milky Way. 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 radiative transfer side, so that the full pipeline becomes numerically demanding even on parallel facilities and then not suitable to make a blind parameters 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 re-calibration involves the two free parameters of GAMESH, namely the star formation efficiency in Ly$\alpha$-cooling halos, $\epsilon_{\ast}$, and the efficiency of supernova-driven winds, $\epsilon_{\rm w}$.

As discussed in LG15, we calibrate the free parameters of the model by requiring the star formation rate, 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 star formation rate 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 modeling strategy adopted to reconstruct the galaxy components (bulge, disk and halo). The interested reader can find in Kennicutt (1998); McKee & Ostriker (2007); Kennicutt & Evans (2012); Bland-Hawthorn & Gerhard (2016) a large collection of critically revised estimates and galaxy modeling techniques.

Besides the different methodologies, in the last years the total stellar mass ($M_{\star})$ inferred for the Milky Way has largely converged to a value of $M_{\star}=[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 colored points in the left panel of Figure 5. In the same panel, we show the redshift evolution of $M_{\star}$, as predicted by the more advanced semi-analytic model described in DB16 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_{\star}=0.09$ and $\epsilon_{\rm w}=0.0016$. The model is in good agreement with the observations, with a final value of $M_{\star}\sim 4.6\times 10^{10}$ M*⊙* for the MW halo candidate. It predicts a total stellar mass of $\sim 6\times 10^{10}$ M*⊙* in the LG and its redshift evolution results consistent with that predicted by DB16 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 DB16 predicts a mass $M_{\star}\sim 7\times 10^{10}$ M*⊙* with a local star formation efficiency $\epsilon_{\star}=0.8$, i.e. one order of magnitude higher than the value required by GAMESH131313Note that in both models the global star-formation efficiency is defined as, $M_{*}/M_{{\rm gas}}=\epsilon_{*}(\Delta t/t_{\rm dyn})$, where $t_{\rm dyn}$ is the redshift dependent halo dynamical time and $\Delta t$ is the time scale assumed for star formation.. The reason for this difference can be ascribed to: (i) the different mass of the final MW dark matter halo, which in DB16 is assumed to be $M_{\rm 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 dark matter evolution of the MW halo predicted by the N-body simulation, with continuous mass exchanges between halos entering the MW merger tree and halos 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\leq 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 DB16.

This can be seen in the middle panel of Figure 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_{\ast}M_{\rm gas}$, while the green dashed line represents to the total $M_{\rm gas}$ predicted by DB16. As illustrated in Section 2.4 and in Figure 2, the early assembly of the MW halo in the N-body simulation is dominated by mass accretion and mergers of nearby halos, and the evolution is similar to the one predicted by DB16 using the semi-analytic merger trees based on the Extended Press Schechter formalism (EPS, Press & Schechter 1974). When $z<4$, the simulated MW-halo grows by many episodes of violent accretion and many minor mergers with halos 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 disk) and a hot halo (coronal) component, as exhaustively detailed in Bland-Hawthorn & Gerhard (2016). The observations reported in the middle panel of Figure 5 refer to the total interstellar medium (ISM) mass in the disk 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_{\rm gas}\sim 1.3\times 10^{11}M_{\odot}$.

Finally, in the right panel of Figure 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 behaviors which reflect their corresponding stellar mass assembly histories, and DB16 predicts a higher metal content at $z=0$, consistent with its higher $M_{\star}(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_{\ast}$ galaxies (Peeples et al., 2014). To be consistent with the observations reported in the other panels of Figure 5, we have computed the mass of metals and dust in the ISM from the fitting functions of Peeples et al. (2014), using a stellar mass in the range $[3-7]\times 10^{10}M_{\odot}$ and we find $M_{\rm Z}=[0.95-4.7]\times 10^{8}M_{\odot}$.

In addition, Peeples et al. (2014) provide an estimate of the CGM metal mass as probed by low- and high-ionization (O vi) species based on COS-halos data and of the mass of dust in the circum-galactic medium (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 gas141414We compute this quantity from equation (24) in Peeples et al. (2014) assuming a stellar mass of $[3-7]\times 10^{10}M_{\odot}$., 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 dark matter halo radial profile at these large radii, we find $M_{\rm Z}(r<150~{}{\rm 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_{\ast}=0.09$ and $\epsilon_{\rm 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 DB16, 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 Figure 6 we show the total star formation rate (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 et al. (1978, SFR$\sim 5$ M*⊙* yr*-1*, dark green cross), McKee & Williams (1997, SFR$\sim 4$ M*⊙* yr*-1*, violet square); Diehl et al. (2006, SFR$\sim 4$ M*⊙* yr*-1*, violet square) to the newest, generally lower values around SFR $\sim 2$ M*⊙* yr*-1* by Robitaille & Whitney (2010, orange minus), Licquia & Newman (2015, light-violet circle), 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 Milky Way (for example 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 errorbars). For the MW-like halo at $z=0$, GAMESH finds a SFR$\sim 4.7$ M*⊙* yr*-1*, a factor of two 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 DB16, and the two SFRs show similar peak values ($\sim 15-20$ M*⊙* yr*-1*), they peak at different redshits. The progressive, quasi parallel decline of the total SFR, results in final SFRs differing by $\sim 1.5$ M*⊙* yr*-1*. As argued for the gas and stellar mass behaviors, we ascribe these discrepancies to intrinsic differences in the merger tree definitions discussed above, MW mass assembly history, particularly at $z<4$, and to the IRA, that naturally introduces an acceleration in the evolution.

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

4 Properties of Milky Way progenitors

So far, we have investigated the global properties of the simulated halos, 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-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 sets based on the FourStar Galaxy Evolution (ZFOURGE) 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 which falls within $\pm 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$) to 6 (in the highest redshift bin, $2.2<z<2.8$).

Figure 7 shows the resulting evolution of the average stellar mass (right panel), SFR (middle panel) and gas fraction (right panel), defined as $f_{\rm gas}=M_{\ast}/(M_{\rm gas}+M_{\ast})$. 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 $\pm 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 which 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\%$ of the MW mass has been built since $z\sim 2.5$. However, the star formation rate and the gas fraction of the simulated galaxies have a shallower evolution in the 3 Gyr period between $z=2.5$ and $z=1$ than found by Papovich et al. (2015). In particular, the peak SFR of $\sim 10\,M_{\odot}/{\rm yr}$ 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 buildup 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 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 panel of Figure 7). The predicted $f_{\rm 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 Instantaneous Recycling Approximation, which affects the efficiency of mechanical and chemical feedback acting on the evolution of individual galaxies, and by the lack of radiative transfer effects. All these different feedback processes, indeed, can strongly affect the evolution of the MW progenitors (e.g. Salvadori et al. 2010b; Salvadori & Ferrara 2012; Graziani et al. 2015).

We have also checked the position of the simulated MW progenitors relative to the galaxy main sequence of star formation. In Figure 8 we show the results using the same redshift bins adopted in Figure 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 Eq. 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_{\ast}<10^{8}M_{\odot}$ show SFRs that can vary by almost one order of magnitude. While the galaxy main sequence can not be constrained by observations in this regime, an increasing scatter towards low stellar masses has already been found in hydro-dynamical 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 $z\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 (Mannucci et al., 2010, FMR) and fundamental plane of metallicity (Hunt et al., 2012; Hunt et al., 2016a, FPZ).

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+\rm{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 star formation rate. 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_{\ast}>10^{9}M_{\odot}$. Similarly to Figure 8, the scatter in the MZR increases with decreasing $M_{\ast}$, a result that appears to be consistent with deep spectroscopic observations which probe galaxies down to $M_{\ast}\sim 3\times 10^{7}M_{\odot}$ at $0.5\leq z\leq 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 mass-metallicity relation proposed by Hunt et al. (2016a, and computed using their Eq. 2 at the average redshift of each bin), shown as the dark grey dashed line, with the shaded region representing a dispersion of $\pm 0.15$ dex.

A similar result is found in Figure 10, where we show the position of the simulated MW progenitors relative to the FMR (left panel) and FPZ (right panel). Here we have reported systems selected in the same redshift bins as in Figure 9. As usual, smaller MW progenitors which populate the lower right 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 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 star formation rate are believed to originate from the interplay between gas accretion, star formation and SN-driven outflows (see, among others, Dayal et al. 2013 and 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 mini- and Ly$\alpha$-cooling halos

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 dark matter halos introduced in Section 2 (see Figure 4) we discuss the properties of the galaxies hosted in mini- and Ly$\alpha$-cooling halos. We first follow their SFR, M*∗, MZ* in time and then we compare their final distribution at $z=0$.

In Figure 11 we show the redshift evolution of their SFR, both in the MW merger tree and in the LG. In the left panel, the total SFR of mini- and Ly$\alpha$-cooling halos 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 halos151515The curves that show the contribution of mini-halos are interrupted at $z=6$ because we assume that SF is suppressed in mini-halos when $z\leq z_{\rm reio}$ (see Section 3).. Despite the number of mini-halos is largely dominant, as shown in Figure 4, their small contribution to the total SFR at all redshifts is due to (i) their low SF efficiency (Salvadori & Ferrara, 2009, 2012; Graziani et al., 2015) compared to Ly$\alpha$-cooling halos, where the SFR is proportional to the halo gas content and controlled by the accretion history (see Figure 3), (ii) their low gravitational potential, which implies an intrinsic inefficiency in accreting gas from the IGM.

Among Ly$\alpha$-cooling halos, a special role is played by the most massive halo which 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 panel of Figure 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 halos 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 onto 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 panel of Figure 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 mini-halos and the sharp drop in their SFR is driven by the increase of the average metallicity of the IGM above $Z_{\rm 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 mini-halos. 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 as the JWST161616 Due to the intrinsic differences in the dynamical configuration of halos in our LG with respect to the observed one, this information should be valuated only in a statistical sense. . 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 et al., 2003; Salvadori et al., 2014; Sarmento et al., 2016). This point will be deeply investigated in the companion work of this paper including full RT and particle tagging (Graziani et al., in prep.).

The remaining redshift evolution is dominated by Pop II stars formed in mini- and Ly$\alpha$-cooling halos (solid thin and thick red lines, respectively). We also show the contribution of Pop II star formation in mini-halos 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_{\rm 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 prep).

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 mass171717Note that the MW galaxy is excluded.. Their distribution follows the galaxy main sequence presented in Figure 8. The three most massive SF galaxies, with $M_{*}>10^{9}M_{\odot}$, have $0.1<$SFR$<1M_{\odot}$ 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}<$ SFR $<10^{-1}$ $M_{\odot}$ yr*-1*, as observed in the Large (LMC) and Small Magellanic Clouds (SMC) which have $M_{\ast}\approx 10^{9}$ $M_{\odot}$, $0.25<$ SFR $<0.63$ $M_{\odot}$ yr*-1* (LMC) and $M\approx 3\times 10^{8}$ $M_{\odot}$, $0.016<$ SFR $<0.039$ $M_{\odot}$ 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$\alpha$-cooling halos with SFR$<10^{-3}$ $M_{\odot}$ yr*-1*, within the detection limits of recent dwarf galaxy surveys, such as the DGS (Madden et al., 2013, 2014).

In the two panels of Figure 12 we show the redshift evolution of the stellar (top panel) and metal (bottom panel) mass , separating the contribution of MW progenitors hosted in mini- and Ly$\alpha$-cooling halos (solid thin and thick blue lines, respectively), and of all the mini-halos found in the LG (dotted black lines). As expected, progenitors in Ly$\alpha$-cooling halos dominate the evolution and the total stellar and metal mass in mini-halos is orders of magnitude smaller, even when all the mini-halos 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$\alpha$-cooling and mini-halos 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 semi-analytic model in which the IRA approximation is relaxed and the inhomogeneous metal enrichment and reionization of the MW environment is accounted for by using analytic prescriptions.

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

Below $z<2$, the oscillating behavior in the thin blue lines are due to continuous mass exchange between small progenitors (both mini and small Ly$\alpha$-cooling halos) orbiting around the MB, before the final merging, which causes the sharp drop at $z=0$181818By definition, the MW merger tree contains only the final MW host halo at $z=0$. A series of mass exchange events through dynamical interactions or destructive mergers involving small Ly$\alpha$-cooling and mini-halos 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 mini-halo population is mainly ascribable to a series of contamination events from highly enriched baryon masses tidally dragged by mini-halos orbiting Ly$\alpha$-cooling halos. An increasing number of small Ly$\alpha$-cooling halos transitioning the mini-halos 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 Milky Way 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 Milky Way-like halo and a significant number of smaller objects and MW satellites. This new simulation allows us to draw the following, main conclusions:

• •

The new DM simulation has the adequate mass resolution to guarantee a reliable Milky Way-type DM halo whose mass and dynamical properties are in agreement with a number of independent simulations and observational estimates. The Local Group surrounding it shows a plethora of intermediate mass Lyman $\alpha$-cooling halos and a vast number of mini-halos, 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.

• •

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 Milky Way observations. We predict a final SFR for the resulting Milky Way system a factor of 2 higher than recent simulations but still in agreement with the many observational uncertainties.

• •

A particle-by-particle reconstruction of the Milky Way 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% of the MW mass has been built since $z\sim 2.5$. However, the star formation rate and the gas fraction of the simulated galaxies have a shallower evolution between $z=2.5$ and $z=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.

• •

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 mass-metallicity relation, fundamental metallicity relation and fundamental plane of metallicity at $0<z<4$ (Mannucci et al., 2010; Hunt et al., 2012; Hunt et al., 2016a). Since these scaling relations are believed to originate from the interplay between gas accretion, star formation and SN-driven outflows, we conclude that the description of these physical processes obtained by GAMESH leads to results consistent with observations.

• •

At all but the highest redshifts, the SFR of the Milky Way is dominated by a multiplicity of galaxies hosted in Lyman $\alpha$-cooling halos, 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 mini-halos in the Local Group is comparable to the SFR along the MW merger tree at $z>6$, indicating that these systems provide an important source of ionizing photons.

• •

Due to efficient metal enrichment, Pop III stars are confined to form in the smallest mini-halos at $z>16$, and their formation rate is larger in the Local Group than along the MW merger tree. This suggests that traces of Pop III star formation are not confined to the MW and its satellites but may be found in external galaxies of the Local Group, although their detection may be challenging even for the next generation of telescopes.

• •

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

• •

The low redshift evolution of all halos, when followed in stellar mass and metal mass, shows the importance of dynamical effects acting onto progenitors which 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 radiative transfer and inhomogeneous metal enrichment with GAMESH.

Acknowledgments

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 & Technology Facilities Council and the Large Facilities Capital Fund of BIS. We also acknowledge PRACE 191919http://www.prace-ri.eu/ for awarding us access to the CEA HPC facility "CURIE@GENCI" 202020http://www-hpc.cea.fr/en/complexe/tgcc-curie.htm 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 n. 306476.

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 halos found in the AQUARIUS, ELVIS, CATERPILLAR and APOSTLE simulations. Table 1 and Table 2 collect halos having M${}_{\texttt{vir}}\sim 1.5-1.9\times 10^{12}$ M*⊙* and summarize their structural properties. The Hamilton/Hall halos 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 satellites distribution, on the other hand, has been extensively studied in the literature and will allow future dedicated comparisons (Mancini et al. in prep).

From Table 2 we infer that our simulation adopts a mass resolution and cosmology compatible with level 13 (Halos Cat-8/LX13 and Cat-12/LX13) in the CATERPILLAR project, while the ELVIS simulation provides halos better resolved by a factor 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 WMAP-7 measurements. Finally, ELVIS has simulated both isolated (iHamilton/iHall) and paired halos (with an M31 analogue at a distance of $\sim 800$ kpc), while the MW halo and the halos identified in the CATERPILLAR sample are isolated.

By checking the values in Table 2 it is evident that the MW halo structural properties, such as its virial radius R${}_{\texttt{vir}}$, concentration parameter $c$ (see Navarro et al. 1996, 2010 for a definition) and maximum circular velocity V${}_{\texttt{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${}_{\texttt{vir}}$, R${}_{\texttt{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 Rs, the size of the "non contaminated region". In our simulation, a cubic volume of 4 cMpc side-length is uncontaminated, which is larger that spherical regions with Rs = 1.1-1.5 cMpc surrounding the other halos. This rich (but still computationally affordable) dynamical information will allow future studies to adopt accurate RT feedback and/or particle-tagging techniques to model in-homogeneous 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 MZ, with independent SAM models coupled to DM simulations. As often discussed in many SAM comparisons projects (Knebe et al., 2015; Hirschmann et al., 2012) 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 halos (Navarro et al., 2010; Terrazas et al., 2016). Despite this intrinsic limitation, when compared with predictions based on Milky Way-like halos 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 panels of Figure 2 and 5 relative to $12<{\rm log}M<13$) or with the FULL case of Figure 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_{\star}=0.02$ on a set of 56 MW-like halos having $3<M_{\star}/10^{10}M_{\odot}<7$ and $0.7<M_{\rm DM}/10^{12}M_{\odot}<3$ (in their jargon $MW-sister$ and $MW-cousins$ halos). They also find SFR$(z=0)=(1-5)$ M${}_{\odot}/yr$. Despite their disk modeling 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 Figure 6 with Figure 9 in their paper).