{'fraction_non_alphanumeric': 0.024702641302502906, 'fraction_numerical': 0.011202360590669922, 'mean_word_length': 4.510215827338129, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Numerical hydrodynamic simulations based on semi-analytic galaxy merger trees: method and Milky-Way like galaxies

We present a new approach to study galaxy evolution in a cosmological context. We combine cosmological merger trees and semi-analytic models of galaxy formation to provide the initial conditions for multi-merger hydrodynamic simulations. In this way we exploit the advantages of merger simulations (high resolution and inclusion of the gas physics) and semi-analytic models (cosmological background and low computational cost), and integrate them to create a novel tool. This approach allows us to study the evolution of various galaxy properties, including the treatment of the hot gaseous halo from which gas cools and accretes onto the central disc, which has been neglected in many previous studies. This method shows several advantages over other methods. As only the particles in the regions of interest are included, the run time is much shorter than in traditional cosmological simulations, leading to greater computational efficiency. Using cosmological simulations, we show that multiple mergers are expected to be more common than sequences of isolated mergers, and therefore studies of galaxy mergers should take this into account. In this pilot study, we present our method and illustrate the results of simulating ten Milky Way-like galaxies since z=1. We find good agreement with observations for the total stellar masses, star formation rates, cold gas fractions and disc scale length parameters. We expect that this novel numerical approach will be very useful for pursuing a number of questions pertaining to the transformation of galaxy internal structure through cosmic time.


INTRODUCTION
Large galaxy redshift surveys, such as the Two Degree Field Galaxy Redshift Survey (Folkes et al. 1999) and the Sloan Digital Sky Survey (Stoughton et al. 2002) have measured and quantified the large scale structure of the Universe. On the other side, computational advances also allow accurate predictions of the distribution of dark matter on large scales based on numerical simulations within the Λ Cold Dark Matter (CDM) paradigm. The predictions of ΛCDM are in excellent agreement with observations on large scales.
On small scales however, it is much more difficult to explain and predict galaxy structures as this requires an understanding of the baryonic processes that shape galaxy moster@mpa-garching.mpg.de properties. The standard picture of hierarchical galaxy formation was introduced by White & Rees (1978) and has been subsequently extended. In this picture, small initial density fluctuations grow with time due to gravitational instability. While the dark matter collapses into haloes with a quasi-equilibrium state through violent relaxation, the baryonic matter falls into the potential wells of these haloes forming a hot gaseous halo for which the self-gravity is balanced by pressure gradients. The gas in this halo is able to cool, which reduces the pressure support and causes the accretion of cold gas from the halo onto a central disc. Here, stars can form depending on the gas density on a characteristic timescale resulting in a rotating stellar disc. A fraction of the newly created stars is short-lived and these stars explode in supernovae (SNe) which can heat the surrounding gas, reducing the efficiency of star formation and in some cases blowing gas out of the galaxy in a galactic wind. Another potential source of feedback is provided by the accretion of gas onto a supermassive black hole (BH) in the centre of the galaxy, called active galactic nuclei (AGN).
Due to the hierarchical nature of the ΛCDM scenario, a dark matter halo constantly accretes new material and eventually other galactic systems merge with it. Such a merger may be accompanied by a strong burst of star formation if the merging galaxies have a similar mass and contained significant amounts of cold gas. At the same time, angular momentum is transferred to the stars in the disc. For a major merger, the orbits of the disc stars are randomised, resulting in the destruction of the discs and the creation of an elliptical galaxy. After such a merger a new gas disc can be created and a new stellar disc formed (Moster et al. 2011). In minor mergers, the disc of the central galaxy typically survives, although it may be thickened and some of the stars can be removed from the disc and transferred to a spherical component while the satellite is completely destroyed (Moster et al. 2010b(Moster et al. , 2012. In this way, galaxy mergers affect the morphology of galaxies and shape their properties. One of the major goals of galaxy formation modeling is to understand the correlations of galaxy internal properties with their formation history and environment. In addition to global properties such as luminosity or stellar mass, star formation rate, and color, we can also study the internal structure and kinematics of galaxies. For example, there are wellknown scaling relations between galaxy luminosity or stellar mass, radial size, and rotation velocity or velocity dispersion (called the Fundamental Plane; the Tully-Fisher relation for late type galaxies or Faber-Jackson for early types are projections of this plane). Moreover, there are strong correlations between galaxy morphological or structural properties (e.g. spheroid vs. disk dominated) and spectrophotometric properties (color or specific star formation rate), which are not yet fully understood.
In order to study the connection between large scale environment and galaxy internal properties, we can specify several requirements for a galaxy formation modeling approach. The first prerequisite is the inclusion of gas and the physical processes to which it is subject. While this may seem obvious, there are many studies in the literature that neglect the gaseous component and use initial conditions that include only the collisionless components of stars and dark matter. However, not only does the gas form stars, subsequently altering the density and gravitational potential, but the gas is also able to radiate energy away through cooling. By this mechanism, the total energy content of a system can be reduced, which is not possible in pure collisionless studies. The next requirement is high enough resolution to capture the internal structure of galaxies. Many cosmological simulations of galaxy formation have ∼ kpc spatial resolution, which clearly does not allow us to say anything about the internal structure or morphology of galaxies. Third, it is important to study the evolution of galaxies within a cosmological background. The cosmological framework has a large impact on the evolution of a galaxy. It specifies when a halo and a galaxy form and how much mass they accrete during their evolution. Furthermore, the cosmology determines the merger histories of the galaxies, i.e. when and how often a galaxy undergoes a merger event. Finally, while it can be useful to study a small number of systems in great detail, in order interpret the significance of any results and make use of the available large observational samples of galaxies, it is helpful to be able to simulate a statistically significant sample of systems, spanning a range of the relevant properties such as halo mass, formation history, etc. In this paper we present a novel approach that allows us to carry out numerical hydrodynamic simulations of galaxies at high resolution, within a proper cosmological context, with much greater computational efficiency than standard methods.

Numerical simulations
Numerical simulations are a powerful tool for studying galaxy formation and evolution. Dissipationless N -body methods, which neglect the gas physics, have been extensively employed. State-of-the-art N -body simulations include the 'Millennium' simulation , the 'Bolshoi' simulation (Klypin et al. 2011), the 'Via Lactea' simulation (Diemand et al. 2007), the 'Aquarius' project (Springel et al. 2008) and the 'GHALO' simulations (Stadel et al. 2009). Collectively, N-body simulations have characterized the formation of dark matter structures from the internal structure of dark matter halos, to the largest structures we expect to form on hundreds of Mpc.
When a dissipational component is included in the simulation, the complexity of the problem increases substantially as does the computational time. One consequence is that cosmological hydrodynamical simulations typically cannot reach the same spatial and mass resolution of N -body runs. Due to the limited resolution, as well as imperfect treatment of feedback processes, cosmological hydrodynamic simulations have long suffered from a set of intertwined problems. They have tended to produce galaxies that are too compact and too bulge-dominated (the angular momentum catastrophe), and with stellar or baryonic masses that are too large compared with their halo masses (the overcooling problem). They do not typically reproduce the red colors and very low star formation rates of observed massive galaxies (another manifestation of the overcooling problem). Recent efforts have demonstrated that with very high resolution, achievable through "zoom" techniques, as well as improved treatment of star formation and stellar feedback, it may be possible to produce disk-dominated galaxies with sizes and baryon fractions that are consistent with observations. (Governato et al. 2004;Robertson et al. 2004;Okamoto et al. 2005;Governato et al. 2010;Agertz et al. 2011;Guedes et al. 2011;Stinson et al. 2012). This is encouraging, but a drawback of this approach is the enormous computational expense. Even with very large expenditures of computational resources, it is only feasible to simulate a handful of systems at this resolution, particularly if one wishes to run the simulations to z = 0.
A third approach neglects the cosmological background and simulates an isolated galaxy or the interaction of two pre-formed galaxies in a binary merger. With this technique, it is possible to study galaxy transformations at very high resolution employing gas physics (Noguchi 1988;Combes et al. 1990;Mihos & Hernquist 1996;Cox et al. 2006;Naab et al. 2006;Robertson et al. 2006;Moster et al. 2011). However, in this approach, the initial conditions (such as the properties of the progenitor galaxies and the orbit) are not derived from a cosmological context, but are specified a pri-ori, based on observations or simply spanning a desired range of values. Perhaps more importantly, on the timescale of the merger process (a few Gyr), one expects a significant amount of new material (gas and dark matter) to be accreted in a cosmological framework, but this accretion has been neglected in most binary merger studies. This accretion can greatly affect the evolution of the galaxies (Moster et al. 2011). Furthermore, many galaxies experience multiple mergers over their lifetime, and indeed, these mergers frequently occur in rapid enough succession that the galaxy does not have time to relax in between (see section 3). Therefore this process may not be accurately simulated via a sequence of binary mergers.

Semi-analytic models
The semi-analytic approach was proposed by White & Frenk (1991), based on ideas presented in White & Rees (1978). This method, which could only predict average quantities, was later extended to model the properties of individual systems, based on merger histories as predicted by the CDM scenario (Kauffmann et al. 1993;Cole et al. 1994;Kauffmann et al. 1994;Baugh et al. 1996;Somerville & Primack 1999;Bower et al. 2006;Croton et al. 2006;Somerville et al. 2008). In a semi-analytic model (SAM) one starts by specifying a cosmological model and then traces the merger histories for a series of dark matter haloes using either an Nbody simulation or the analytic extended Press-Schechter (EPS) formalism. The evolution of the baryonic component in the haloes is followed by using simple, yet physical analytic recipes for gas cooling and accretion, star formation and feedback. Cooling converts hot gas into cold gas, star formation converts cold gas into stars, and feedback converts cold gas into hot gas and in some cases ejects it from the halo. When two haloes merge, it is assumed that the galaxies merge on a time scale set by dynamical friction, possibly altering the properties of both systems. A sample of model galaxies that represents the galaxy population in the present-day Universe can then be built by modelling a large number of haloes that sample the halo mass function. The resulting star formation and chemical evolution histories can be convolved with stellar population synthesis models to predict observables such as luminosity and colour. For an excellent recent review on SAMs we refer to Baugh (2006).
Thus, in a SAM the complicated astrophysical processes that are responsible for the formation and evolution of galaxies are modelled as a set of recipes which carry a number of free parameters (it should be kept in mind, however, that numerical hydrodynamic simulations carry the same free parameters to characterize the sub-grid physics). As the current understanding of many of these processes is limited, the model parameters cannot be derived from first principles. Instead, the model is normalised with a set of observational constraints, and with the help of more detailed hydrodynamical simulations. The advantage of SAMs is their flexibility and their computational efficiency. It is possible to run many realisations in a short amount of time, and thereby test the effects of the various assumptions and model parameters. It is also possible to create 'mock catalogs' of large samples of galaxies (millions), which can be compared with large observational samples. Moreover, one can include more physi-cal processes, such as AGN feedback, albeit in a schematic way. Currently, SAMs are much more successful than numerical cosmological simulations at reproducing the global properties of galaxies. Furthermore, SAMs reproduce quite well the results of hydrodynamic "zoom" simulations when similar physics is included (Hirschmann et al. 2012).
The main drawback of SAMs is that they yield only approximate and schematic information about the internal structure and morphology of galaxies. For example, most SAMs provide an estimate of the ratio of the spheroid mass to the disc mass, based on the merger history of the galaxy and the assumption that near-equal mass mergers destroy discs and build spheroids (e.g. Fontanot et al. 2011). Many SAMs also include a simple recipe for computing the radial size and rotation velocity of the disc component based on angular momentum considerations (e.g. Somerville et al. 2008), and recent work provides an approach for predicting the size and velocity dispersion of spheroids formed in mergers (Covington et al. 2011). However, in the existing SAMs, these recipes are largely based on binary merger simulations, which we have argued may not accurately capture the associated physics because of the neglect of the cosmological accretion. In summary, we may consider the SAM to be an efficient tool for predicting the global properties of large samples of galaxies, but of limited use for studying galaxy internal properties.

Combining simulations and semi-analytic models
In this paper, we present a new approach that leverages the complementary strengths of semi-analytic models and hydrodynamic merger simulations to efficiently simulate the evolution of galaxies at high resolution and within a cosmological context. We start with a dark matter halo merger tree extracted from a dissipationless N-body simulation. We run the SAM within this halo merger tree to construct a galaxy merger tree, which specifies when each galaxy enters a larger halo and becomes a satellite, and the global properties of this galaxy (e.g. stellar mass of a bulge and disc component, cold gas mass, radial size) when it enters the halo. We use these SAM predictions to specify the initial conditions for a sequence of mergers, which we simulate with a full numerical hydrodynamic code. Specifically, we evolve the galaxy in the main branch with our hydrodynamical code, and include satellite galaxies in the simulation at the time when they cross the virial radius of the larger halo (with this time specified by the Nbody similulation) using the satellite properties predicted by the SAM. In addition, we include the growth of the main halo due to cosmological accretion, as specified by the merger tree. This approach is particularly well suited for studies that focus on the evolution of the internal structure of galaxies, including different components such as a thin and thick disc, spheroid, and stellar halo. While many applications of our model focus on the evolution of the central galaxy, we can also address the evolution of the satellite population. Finally, it is our goal to develop generalised semi-analytic prescriptions based on the results of our simulations. In this way, the existing semi-analytic recipes can also be improved.

METHODS
In this section we briefly describe the tools that are used in this paper. These are the simulation code GADGET-2 that was employed both for the cosmological N -body simulation from which our merger trees were extracted and the hydrodynamic merger simulations, the code to create initial galaxy models, and the semi-analytic model used to populate N-body merger trees with galaxies. Throughout this paper we adopt cosmological parameters chosen to be consistent with results from WMAP-3 (Spergel et al. 2007) for a flat ΛCDM cosmological model: Ωm = 0.26, ΩΛ = 0.74, h = H0/(100 km s −1 Mpc −1 ) = 0.72, σ8 = 0.77 and n = 0.95. These parameters differ only slightly from the currently favored WMAP7 parameters (Komatsu et al. 2011). We adopt a Kroupa (2001) IMF and compute all stellar masses accordingly.

N-body and Smoothed Particle
Hydrodynamics code For all simulations in this work we employ the parallel TreeSPH-code GADGET-2   Monaghan 1992) to evolve the gas using an entropy conserving scheme (Springel & Hernquist 2002). Radiative cooling is implemented for a primordial mixture of hydrogen and helium following Katz et al. (1996), and a spatially uniform time-independent local UV background in the optically thin limit (Haardt & Madau 1996) is included. All simulations have been performed with a high force accuracy of α force = 0.005 and a time integration accuracy of ηacc = 0.02 (for further details see . Star formation and the associated heating by supernovae (SN) is modelled following the sub-resolution multiphase ISM model described in Springel & Hernquist (2003). The ISM in the model is treated as a two-phase medium with cold clouds embedded in a hot component at pressure equilibrium. Cold clouds form stars in dense (ρ > ρ th ) regions on a timescale chosen to match observations (Kennicutt 1998). The threshold density ρ th is determined self-consistently by demanding that the equation of state (EOS) is continuous at the onset of star formation. SN-driven galactic winds are included as proposed by Springel & Hernquist (2003). In this model the mass-loss rate carried by the wind is proportional to the star formation rate (SFR)Ṁw = ηṀ * , where the mass-loading-factor η quantifies the wind efficiency. Furthermore, the wind is assumed to carry a fixed fraction of the supernova energy, such that there is a constant initial wind speed vw (energy-driven wind). We do not include feedback from accreting black holes (AGN feedback) in our simulations.
We compute the parameters for the multiphase feedback model following the procedure outlined in Springel & Hernquist (2003) in order to match the Kennicutt Law. For a Kroupa IMF the mass fraction of massive stars is β = 0.16 resulting in a cloud evaporation parameter of A0 = 1250 and a SN "temperature" of TSN = 1.25 × 10 8 K. Finally, the star formation timescale is set to t 0 * = 3.5 Gyr. For the galactic winds we adopt a mass-loading factor of η = 1 and a wind speed of vw ∼ 500 km s −1 .

Merger Trees
For this study we use dark matter merger trees drawn from an N -body simulation run with the GADGET-2 code. The initial conditions for the WMAP-3 cosmology were generated using the GRAFIC software package (Bertschinger 2001). The simulation was done in a periodic box with a side length of 100 Mpc, and contains 512 3 particles with a particle mass of 2.8 × 10 8 M and a comoving force softening of 3.5 kpc. Starting at a redshift of z = 43, 94 snapshots were stored until z = 0, equally spaced in expansion factor (∆a = 0.01).
Dark matter haloes were identified in the simulation snapshots using a Friends of Friends (FOF) halo finder with a linking parameter of b = 0.2. Substructures inside the FOF groups are then identified using the SUBFIND code. For the most massive subgroup in a FOF group the virial radius and mass are determined with a spherical overdensity criterion using the fitting function by Bryan & Norman (1998). The minimum particle number for haloes is set to 20, resulting in a minimum halo mass of 5.6×10 9 M . The number of parent haloes found at z = 0 with this method is ∼ 135 000 and parent halo masses range between 10 10 M and 10 15 M . The merger trees are constructed for all parent haloes at z = 0 by connecting haloes between the 94 catalogues. The branches of the trees were determined by linking every halo to its most massive progenitor at previous snapshots. In total, we have 41 000 merger trees.

Galaxy models
To construct the galaxy models used in our simulations we apply the method described in Springel et al. (2005) with the extension by Moster et al. (2011). Each system is composed of a cold gaseous disc, a stellar disc and a stellar bulge with masses Mcg, M disc and M b embedded in a halo that consists of hot gas and dark matter with masses M hg and M dm .
The gaseous and stellar discs are rotationally supported and have exponential surface density profiles. The scale length of the gaseous disc rg is related to that of the stellar disc r d by rg = χr d . The vertical structure of the stellar disc is described by a radially independent sech 2 profile with a scale height z0, and the vertical velocity dispersion is set equal to the radial velocity dispersion. The gas temperature is fixed by the EOS, rather than the velocity dispersion. The vertical structure of the gaseous disc is computed self-consistently as a function of the surface density by requiring a balance of the galactic potential and the pressure given by the EOS. The spherical stellar bulge is non-rotating and is constructed using the Hernquist (1990) profile with a scale length r b . The dark matter halo has a Hernquist (1990) profile with a scale length rs, a concentration parameter c = rvir/rs and a dimensionless spin parameter λ.
The hot gas is modelled as a slowly rotating halo with a spherical density profile following Moster et al. (2011). We employ the observationally motivated β-profile (Cavaliere & Fusco-Femiano 1976;Jones & Forman 1984;Eke et al. 1998): which has three free parameters: the central density ρ0, the core radius rc and the outer slope parameter β. We adopt β = 2/3 (Jones & Forman 1984), rc = 0.22rs (Makino et al. 1998) and fix ρ0 such that the hot gas mass within the virial radius is M hg . The temperature profile is fixed by assuming an isotropic model and hydrostatic equilibrium inside the galactic potential such that it is supported by pressure. In addition, the hot gaseous halo is rotating around the spin axis of the disc with a specific angular momentum j hg , which is a multiple of the specific angular momentum of the dark matter halo j dm such that j hg = αj dm . The angular momentum distribution is assumed to scale with the product of the cylindrical distance from the spin axis R and the circular velocity at this distance: j(R) ∝ R vcirc(R).
High resolution cosmological simulations that succeed in reproducing disc sizes find that at low redshift (z ∼ < 2), α is greater than unity. This is because feedback processes preferentially remove low angular momentum material from the halo (Governato et al. 2010). Similarly, Moster et al. (2011) constrain α by using isolated simulations of a MWlike galaxy and demanding that the evolution of the average stellar mass and scalelength found observationally be reproduced. The model that agrees best with the observational constraints has a spin factor of α = 4. We use this value throughout this work.

Semi-analytic model
In this paper we use the SAM constructed by Somerville & Primack (1999) in its current implementation (Somerville et al. 2008. In this model each dark matter halo is assigned two properties: the spin parameter λ and the concentration parameter c. Each top-level halo is assigned a value of λ by selecting values randomly from a lognormal distribution with meanλ = 0.05 and width σ λ = 0.5 (Antonuccio-Delogu et al. 2010). The initial density profile of each halo is described by the NFW form. For N -body merger trees the concentration parameter and the halo position are extracted from the simulation. After accretion the orbital evolution of each subhalo is computed using the dynamical friction formula given by Boylan-Kolchin et al. (2008), including the effects of tidal stripping and destruction.
The SAM adopts a simple but fairly standard spherical cooling model, which makes use of the multi-metallicity cooling function of Sutherland & Dopita (1993). The hot gas density profile is assumed to be that of a singular isothermal sphere. In order to model the effects of a photo-ionising background, the SAM follows Gnedin (2000) and Kravtsov et al. (2004) and defines a 'filtering mass' MF which is a function of redshift and depends on the re-ionisation history of the Universe. The cold gas profile is assumed to be an exponential disc with a scalelength proportional to the scalelength of the stellar disc. Only gas lying above a critical surface density threshold is available for star formation. In the 'quiescent' phase star formation is based on the empirical Schmidt-Kennicutt law with the appropriate normalisation for a Kroupa IMF. Massive stars and supernovae impart thermal and kinetic energy to the cold interstellar medium, and cold gas may be "reheated" from the disc and either deposited in the hot halo or ejected from the halo altogether.
When a subhalo loses its orbital energy due to dynamical friction the satellite merges with the central galaxy, driv-ing a starburst. The efficiency of star formation in a mergertriggered burst is parametrized as a function of the mass ratio of the merging pairs and the gas fraction of the progenitors and calculated with the model proposed by Hopkins et al. (2009), which is based on binary merger simulations. It is assumed that during every merger with a mass ratio above µ = 0.1, a fraction of the disc stars is transferred to the spheroidal or bulge component, again following the prescription outlined in Hopkins et al. (2009).
The SAM contains recipes for black hole growth and AGN feedback due to both "bright mode" AGN-driven winds and "radio mode" heating by radio jets (Somerville et al. 2008). However, at the halo mass scale considered here, AGN feedback has little impact on the predictions.

MULTIPLE MERGERS IN A ΛCDM UNIVERSE
In order to help motivate this work we first address whether most mergers can be treated as "isolated" or whether multiple mergers are expected to be cosmologically significant.
In an 'isolated' merger (also called 'binary' merger), two galaxies merge and become dynamically relaxed before the remnant merges with another galaxy. In 'multiple' mergers, two galaxies have not had enough time to merge and become relaxed when another galaxy already enters their common halo. If isolated mergers are the standard event, then it is possible to study single merger events with specified parameters (e.g. in simulations) and then 'stack' these merger events for a given merger history, in order to determine the effects of the mergers on the galaxy. This can be seen as a linear process. However, if multiple mergers are more common, then it may not be possible to string the results from binary merger simulations together. One then has to study the events with three or more galaxies involved which will have a different impact on the galaxy properties than two or more binary mergers. The merger process then becomes non-linear and the (already large) parameter space for the merger parameters multiplies. In order to study whether binary or multiple mergers are more common, we make use of the cosmological merger trees described in Section 2.2. We divide the merger trees according to the mass of their parent halo at z = 0 with a bin size of ∆ log(M/M ) = 0.5. We then identify all merger trees of these systems with two or more halo mergers and for every pair of mergers we record the time difference between the halo mergers. If e.g. the halo of the satellite galaxy A merges with the halo of the central galaxy C at a cosmic time t1 = 8.0 Gyr and later the halo of the satellite galaxy B merges with the halo of the central galaxy at t2 = 10.0 Gyr, we store the time difference ∆t = 2.0 Gyr. We divide this time by the dynamical time of the halo t dyn = rvir/Vvir at the time of the first halo merger t1 and distribute the resulting number in bins of width ∆t/t dyn = 1. Finally, we count the number of mergers pairs with a given mass ratio sequence and since a given redshift for each time bin. The sequence 1:4 → 1:10 for example means that the first merger has a mass ratio above 1:4 and the second merger has a mass ratio above 1:10. In order to get the probability for a merger sequence we normalise this number by dividing it by the number of mergers with the mass ratio of the first merger. We show the results of this analysis in Figure 1. The probability for two mergers happening within a few dynamical times of the halo is always higher than the probability for a larger time difference. This indicates that usually the second satellite galaxy enters the parent halo before the first satellite had time to merge with the central galaxy, which happens only after several dynamical times, depending on the mass ratio of the merger. For the mergers that happen since z = 1 (upper row in the figure) his trend is most obvious. Here, the second satellite enters the halo within three dynamical times after the first satellite entered for more than 50 per cent of all merger pairs independent of mass ratio. If we start at a higher redshift, the merger pairs are on average a little more separated. Still for more the 50 per cent of all merger pairs, the second satellite enters the halo within five dynamical times after the first one since z = 5.
In summary, we find that multiple mergers are more common than sequences of isolated binary mergers. Therefore it is not optimal to focus on binary mergers and afterwards try to string their effects together. Thus merger simulations which consider several satellites entering the halo and orbiting together within the halo have to be performed and analysed. If parameters for the galaxies and the orbits are chosen from a multidimensional grid, the parameter space is then too large to cover. This means that one has to reduce this huge parameter space by selecting only those parameters that are common for the chosen cosmology. N -body merger trees populated with galaxies using a SAM offer this option and are therefore very practical to use as initial conditions for the merger simulations. In this way, mergers that are common in the Universe are automatically selected. Thus, the important regions of the parameter space are naturally covered.

SIMULATIONS OF GALAXY MERGER TREES
Galaxy merger trees obtained from N -body merger trees in combination with a SAM provide natural initial conditions for merger simulations with high resolution. However, if we compare the galaxy models predicted by the SAM to the model galaxies used in merger simulations, we immediately recognise an important difference: the standard approach in binary merger simulations neglects the accretion of material coming from outside the initial virial radius. Standard codes used to create model galaxies only employ a dark matter halo of fixed mass. This means that all material that falls into the halo which is not included in the merger tree (i.e. does not come through satellites) is not accounted for.
On the other hand, cosmological N-body simulations show a substantial increase of the mass of galactic dark matter haloes through "smooth accretion", that then must be taken into account in a cosmologically based study of galaxy evolution. We now describe our method for including this smooth accretion in our simulations.

Modelling the growth of the dark halo
In order to construct a merger tree one must define a lower mass limit (mass resolution) below which the tree is truncated. This means that all haloes that are smaller than this mass limit and fall into a larger halo are not considered as merging haloes, but as so-called 'smooth accretion'. Thus, even if there are no mergers in a given time-step, the parent halo gains mass. In (non-cosmological) simulations of galaxies or mergers this smooth accretion is usually neglected.
We model this accretion by placing the additional dark matter particles into small spherical systems which will be denoted as 'DM-spheres' in the following. We place all DMspheres at a specific distance from the halo centre such that the accretion rate of the DM-spheres matches the smooth accretion rate of the merger tree. The DM-spheres are uniformly distributed around the halo. In order to prevent the DM-spheres from falling to the very centre of the halo and thereby disturbing the disc, we give the DM-spheres an initial velocity in a random direction orthogonal to the radius vector. This of course, reduces the initial distance to the centre compared to freely falling particles. As we require a density profile that leads to a quick dissolving of the DMsphere once it enters the halo (since the DM-spheres are not meant to represent substructure but smooth accretion) we employ a three dimensional Gaussian density profile. Note, that we do not add gas particles to the DM-spheres.

Creating a spherical particle distribution
We model each DM-sphere with a three dimensional Gaussian density profile where MDMS is the total mass of the DM-sphere and rDMS is a scalelength. The mass M (r) enclosed within a radius r is then found by integrating the density profile and the gravitational potential is found by integrating Poisson's equation.
Using N particles of mass mi for each DM-sphere, a Gaussian density profile can by simply achieved by drawing a random number from a Gaussian distribution with width rDMS/ √ 3 for each cartesian coordinate. Deriving the velocity for each particle is more complicated and has to be done using the distribution function as a function of the relative energy f (E) which can be obtained with the so-called 'Eddington inversion' (see e.g. Binney & Tremaine 1987).
In order to solve for f (E), we create a logarithmically spaced array in radius with 10 5 bins. We fix the minimum and maximum of the bins at 10 −3 and 10 2 times the scalelength. On this finely spaced grid, we define the values for ρ(r), M (r) and Ψ(r) and obtain the derivatives by finite differencing. The distribution function is then obtained by numerically integrating the Eddington equation. We find that it is well approximated by the fitting function where Q = E/Ψ0 and Ψ0 = −Φ0 = 2/πGMDMS/rDMS is the maximum relative potential. With the distribution function, we can now find the velocity for each particle with a rejection sampling technique. The direction of the velocity is randomly chosen from the unit sphere. We test this model by creating N -body realisations of a DM-sphere of mass MDMS = 10 9 M . This mass is below the resolution limit of the cosmological simulation used in section 3 and is thus accounted for as smooth accretion. We model the DM-sphere with N = 100, 200 and 1000 particles and a scalelength of rDMS = 2, 5 and 10 kpc. We let each realisation evolve in isolation for 5 Gyr and measure the scalelength at every time-step. The results of this analysis are shown in Figure 2. All models are stable over the 5 Gyr of evolution, independent of the initial scalelength. For N = 100 the scalelength increases by about 50 per cent until the end of the simulation, while for larger particle numbers the scalelength deviates only slightly from the initial value. For all particle numbers, models with a larger initial scalelength are more stable.

Placing the dark matter systems around the halo
With the model presented in the last section we can create DM-spheres that are stable in isolation, but due to their shallow potential they dissolve quickly when orbiting a dark matter halo. The next step in modeling the smooth accretion of a dark matter halo is to place dark matter at the right position, such that the accretion history of the halo is reproduced. Thus, the main difficulty here is where to position the DM-spheres and with which mass.
As we want to simulate merger trees, we first extract the mass accretion history of the 'main branch' from the Nbody simulation, i.e. the mass of the main halo and its most massive progenitors as a function of cosmic time. The next step is to choose a starting redshift zstart which corresponds to a starting time tstart, at which we start the simulation. If we assume that the profile of the parent halo does not change, which is true for our initial conditions, the virial mass of the parent halo increases as the background density decreases towards lower redshift (for a lower ρcrit the radius that contains a certain overdensity becomes larger). We thus subtract the increasing virial radius at every time-step such that we are left with the accreted mass. Furthermore, we subtract the mass that has been accreted through mergers of resolved smaller haloes, as those will be explicitly simulated as merging haloes. We are thus left with the mass M smooth (t) that has been added to the parent halo due to smooth accretion.
In the next step we choose the number of DM-spheres NDMS, which we set equal to the number of time-steps we will use for the model. We then interpolate M smooth (t) for every time-step ti using a Bezier curve which traces the mass accretion history and is equal to the mass found in the merger tree at tstart and tN DMS . We then compute the mass that is accreted within every time-step MDMS(ti), and use this mass for the DM-sphere that corresponds to this time-step. If this mass is smaller than zero (i.e. the parent halo loses mass), we set the mass of the corresponding DMsphere to zero (i.e. the DM-sphere is omitted). However, we store this lost mass and set the mass of the DM-spheres that correspond to the following time-steps to zero, until the lost mass is balanced by accreted mass. Thus the overall mass growth via smooth accretion over the whole accretion history is reproduced. Finally, we set the particle mass for each DM-sphere particle equal to the mass of the dark matter halo particles m dm and compute the number of particles of each DM-sphere: Ni = MDMS(ti)/m dm .
Having determined the mass of each DM-sphere, we now need to specify where to place them. For this we make use of the free-fall time for a point mass at a radius r: t ff = π[r 3 /8GM (r)] 1/2 , where M (r) is the total mass enclosed within r. The time it takes a DM-sphere to enter the virial radius of the halo tenter = ti − tstart is shorter than the free-fall time and only slightly larger than the free-fall time at radius r minus the free-fall time at the virial radius rvir, as the DM-sphere already has a velocity towards the centre when crossing the virial radius. We can thus assume that tenter has a value of tenter ∼ > t ff (ri) − t ff (rvir), where ri is the initial distance of the DM-sphere from the halo centre and rvir is the virial radius of the halo at ti. The initial distance of each DM-sphere is then given by and we place the centre of each DM-sphere randomly on a sphere with radius ri. In order to prevent the DM-spheres from passing through the centre of the halo, we assign an initial velocity to every DM-sphere in a random direction orthogonal to the radius vector. This leads to a larger time taken for a DM-sphere to enter the virial radius which partly compensates for smaller estimate of tenter.
We are thus left with three free parameters: the number of DM-spheres NDMS, their scalelength rDMS, and their initial velocity vinit. In order to fix these values we select a merger tree from the simulation box, presented in section 3, that has no mergers above a mass ratio µ = 0.02. We can thus assume that all mass accreted by the parent halo is from smooth accretion. The index number of this merger tree is 1811 and the virial mass at z = 0 is 10 12 M . We start the simulation at zstart = 1, where the parent halo has a virial mass of 6.5 × 10 11 M , a virial velocity of Vvir = 150 km s −1 , We construct an N -body realisation of this dark matter halo with the initial conditions generator as described in section 2.3, but omitting all baryonic components. We employ N halo = 200 000 particles, and add dark matter particles in DM-spheres to this model with the method outlined above. For this we choose the fiducial parameters as NDMS = 400, rDMS = 5 kpc and vinit = 30 km s −1 , and vary each parameter around these values, while keeping the other two parameters fixed. For about 400 DM-spheres the number of particles per DM-sphere is ≈ 200 which we have shown to be sufficient for stable DM-spheres in isolation. We evolve each system until z = 0 with a softening length of = 0.7 kpc.
For each simulation we measure the virial mass as a function of cosmic time using a spherical overdensity criterion and the fitting function by Bryan & Norman (1998). The results are presented in Figure 3, where the virial mass measured in the simulations is compared to that found in the merger tree. For our fiducial parameters, we are able to reproduce the accretion history of the halo very well. In the left panel the number of DM-spheres is varied. We see that neither a lower nor a higher value than the fiducial one results in a difference for the mass that can be accreted. In the middle panel we vary the scalelength of the Gaussian profile. Also for this parameter we see that the result is independent of its value. In the right panel the initial velocity is varied. Values of vinit = 0 and 30 km s −1 lead to a very similar mass accretion history. For vinit = 60 km s −1 the DM-spheres enter the virial radius slightly later, as the orbital energy is higher and the DM-spheres need to lose some angular momentum first. For higher values (vinit ∼ > 90 km s −1 ) the virial mass of the halo is too low, compared to the merger tree. The reason is that for some DM-spheres the initial velocity is high enough, such that they orbit within the virial radius  of the parent halo just for a short time, or even not at all. Instead their pericentric distance has the same dimension as the virial radius. Thus they do not fall towards the centre, but orbit the halo on a large distance where the density is not high enough such that they can lose their angular momentum due to dynamical friction. This means that the initial velocities of the DM-spheres should not be larger than half of the virial velocity of the halo. Altogether we find that our model is able to reproduce the smooth accretion history of a dark matter halo very well and does not depend on the exact values of the model parameters.
We finally need to check whether the Gaussian DMspheres remain bound when they enter the parent halo, or if they are quickly dissolved such that they can be seen as smooth accretion as intended. For this we compute surface density maps of the system at six redshifts. In Figure 4 the maps are plotted on the scale of 1 Mpc for the model with the fiducial parameters. As the virial radius is of the order of 200 kpc (depending on redshift), the total halo is shown. In the first panel at z = 1, all 400 DM-spheres are clearly visible, as their density is larger than the surrounding background density. As time elapses, the DM-spheres fall towards the centre of the halo. When the DM-spheres enter the halo, they are able to stay bound for a short amount of time, corresponding to less than half an orbital period, and are then dissolved. In the last panel at z = 0, most of the DMspheres have been destroyed. Only those DM-spheres that have just entered the halo are still bound as they have not been within the halo long enough to be tidally destroyed.
As most of the DM-spheres can remain bound for almost SAM Ci S1 S2 S3 S4 Cf ICs Simulation dark matter hot gas cold gas stellar disc stellar bulge ti t1 t2 t3 tf time Figure 6. Schematic view of the combination of semi-analytic models and merger simulations. On the left side the semi-analytic merger tree is shown, with time running from top to bottom. At the starting time t i , a numerical realization of the central galaxy is created using the properties of the central galaxy in the SAM at this time. This system is then simulated with the hydrodynamical code until the time of the first merger t 1 , where the first satellite S 1 galaxy is put into the simulation using the properties predicted by the SAM. The resulting merger is simulated until the next galaxies (S 2 and S 3 ) enter the main halo at t 2 , at which point they are also included in the simulation. This procedure is repeated for all mergers above a specified mass ratio until the final time of the run t f .
half an orbital period, they are able to reach the pericentre of their orbit before they are dissolved. If the DM-spheres are freely falling, this would be problematic for a central disc of stars, since the constant bombardment of small objects can disturb the disc and possibly thicken it. Therefore we modelled the blobs with an initial velocity, such that the pericentre is large enough and the DM-spheres do not perturb the disc. In Figure 5 the surface density maps are plotted on a scale of 50 kpc such that the region within half of the Hernquist scalelength of the halo is shown. This is the typical scale up to which exponential discs with scalelengths of a few kpc can be resolved. As we can see, the density profile remains the same throughout the simulation, and no DMspheres can be identified in this region as here the density of a DM-sphere is lower than the density of the halo. This shows that our method, which models the smooth accretion with small bound objects instead of a uniform distribution of dark matter particles, does not lead to a perturbation of central objects.

Simulations of Semi-Analytic Merger Trees
In this section we describe how we use the SAMs to populate N -body merger trees with galaxies and then use them as the initial conditions for hydrodynamical multiple merger simulations. A schematic view of our method is presented in Figure  6. In a first step we select a dark matter merger tree from the large-scale N -body simulation and use the SAM to predict the properties of the baryonic components of each halo at every timestep. A resulting galaxy merger tree is shown in the left side of Figure 6 with time running from top to bottom. We then choose a starting time ti from which we want to simulate this tree. In our simple example, the main system experiences four mergers after ti. We then use the predictions of the SAM for the central galaxy of the main halo at ti and create a particle realisation with the galaxy generator, as indicated by the brown arrow. This model galaxy is shown in the top right of Figure 6 and in our example consists of a dark matter halo (grey), a hot gaseous halo (red), a cold gaseous disc (blue), a stellar disc (yellow) and a small stellar bulge (green). We evolve this galaxy with our hydrodynamical code until the time t1 when the first satellite galaxy S1 enters the main halo. A particle realisation of the S1 system (dark halo and galaxy) is then created using the semi-analytic prediction for the galaxy properties and included in the simulation at the virial radius of the main halo. The orbital parameters of S1 (position and velocity at the time of accretion) are taken directly from the N -body simulation. We evolve this merger with the hydrodynamical code, until the next satellite galaxy enters the main halo. In our example, two galaxies enter the halo at the same time, and thus we create particle realisations using the predictions of the SAM for both galaxies S2 and S3 and include them in the simulation by positioning them at the virial radius. This procedure is repeated for all merging satellites until the final time t f . In this way, we naturally include multiple mergers, when an already merging galaxy has not been fully accreted as the next galaxy is entering the halo. This is indicated at t3, where the galaxy S3 is still orbiting while the satellite S4 enters the halo. In the example shown in Figure 6, the central galaxy grows both through accretion of dark matter and gas (e.g. from ti to t1 or from t3 to t f ) and through mergers of satellite galaxies (e.g. from t1 to t2).
Our approach thus requires two main steps: creating particle realisations of galaxies as predicted by the SAM and combining these initial conditions in simulations as determined by the merger tree (i.e. every galaxy has to enter the simulation at the specified time and position). This means we first have to specify how the information about the galaxy properties that is computed with the SAM is transformed into three dimensional particle based galaxy models that can be simulated with a hydrodynamic code. Then we have to determine how the satellite galaxies are included in the simulation, so their position and velocity have to be calculated.
We note that the starting time ti, or the starting redshift zi, respectively, can also be chosen such that the simulation starts at a very early epoch. In this case the central galaxy would consist only of a dark matter halo and hot gas in this halo. The dark matter will then grow by mergers and smooth accretion, and the stellar disc and bulge will form as a result of cooling and accretion of gas, and merger events.

Creating particle realisations of semi-analytic galaxies
In order to create particle realisations of the semi-analytic galaxies, we need to specify those galaxies that are included in our simulation. For this we first have to choose a starting redshift zi and a minimum merger ratio µmin. We define this ratio as the mass of the dark matter in the entering subhalo divided by the dark matter mass of the main halo at the time the satellite passes the virial radius. The starting redshift and the minimum merger ratio are in principle free parameters and can be chosen according to the physical problem one wants to address. This makes our approach very flexible.
In a semi-analytic merger tree we identify the central galaxy at the starting redshift, all satellite galaxies within the virial radius at this time (which have not merged with the central galaxy yet) and all satellite galaxies that enter the main halo at a later time. From these galaxies we only select those that fulfill our merger mass ratio criterion. For every selected galaxy we record the value for the time the galaxy enters the main halo tenter, the virial mass of the dark matter halo Mvir, its concentration c, its spin parameter λ, and the semi-analytic predictions for the masses of the hot gaseous halo M hg , the cold gaseous disc Mcg, the stellar disc M disc, * and the stellar bulge M bulge , and the scalelength of the stellar disc r disc . All these quantities are taken at ti for the central galaxy and satellites that are already with the main halo at ti, and at tenter for all satellite galaxies that are entering the main halo after ti. Additionally, for all satellite galaxies we record the dynamical friction time.
The remaining structural parameters of every galaxy that are not directly predicted by the SAM are based on empirical scalings. The scaleheight of the stellar disc is assumed to be a fraction of its scalelength z0 = ζr disc , where typically ζ = 0.15. Similarly the scalelength of the gaseous disc is related to that of the stellar disc by rgas = χr disc , with a standard value of χ = 1.5. The core radius rc of the hot gaseous halo is related to the scalelength of the dark matter halo rs by rc = ξrs with a fiducial value of ξ = 0.22 and the slope parameter of the β-profile is set to β hg = 2/3.
In order to select the number of particles in every component of the galaxies, we determine the semi-analytic prediction for the final stellar mass of the central galaxy M * ,f and choose a number of stellar particles in the central galaxy N * that we would like to obtain at the end of our simulation in order to get the stellar particle mass m * = M * ,f /N * . As the simulation code produces Ng generation of stellar particles from every gas particle due to star formation (where Ng is typically 2), we set the mass of every gas particle to mgas = Ngm * . In this way, all stellar particles (old and new) have the same mass. The mass of the dark matter particles is finally selected as m dm = κm * where κ is a free parameter. If κ is chosen too low, the dark matter particles have a very low mass which results in a very large number of dark matter particles and thus in a high computational cost. For high values of κ the mass of the dark matter particles can become too large, such that these massive particles perturb the disc component which results in numerical disc heating. Typically values of κ ∼ 15 are both computationally efficient and lead to no measurable heating. In order to prevent a component from being unstable due to a low number of particles, we remove a component if its number of particles is lower than a minimum particle threshold Nmin.
Additionally, for simulations where we are interested in the satellite population and not in the central galaxy, we can allow for a higher resolution in the satellite galaxies. This is done by dividing the mass of the particles of every satellite by a number Nres,sat. This number, however, should not be chosen to be too high, since the particle masses of the satellite would be much lower than the particle masses of the main system. This can lead to mass segregation, which means that the particles of higher mass (i.e. the central galaxy particles) will preferentially settle at the bottom of the potential. We wish to avoid this numerical effect which requires that Nres,sat ∼ < 10. The last parameter that needs to be fixed is the gravitational softening length . In order to ensure that the maximum gravitational force exerted from a particle is independent of its mass, we scale the softening lengths of all particle species (dark matter, gas and stars) with the square root of the particle mass (Dehnen 2001). The normalisation of this relation is obtained with a free parameter 1 that specifies the softening length for a particle of one internal mass unit of the code (i.e. 10 10 M ). The softening length is thus given as = 1 mpart/10 10 M , where mpart is the mass of the particle. For our fiducial choice of 1 = 32 kpc and a typical particle mass of 10 5 M this results in a softening length of 100 pc.
Using this recipe, we create particle realisations of all selected galaxies. These can now be used as initial conditions for the hydrodynamical simulation. The next task is thus to specify when a galaxy enters the main halo, and at which positions and with what velocity.

Performing the multiple merger simulation
We start with the particle realisation of the central galaxy and move into its rest frame. For all satellite galaxies we have to compute the relative position and velocity with respect to the central galaxy at the time when they are included in the simulation. Since we use merger trees drawn from simulations, we can directly extract the relative initial positions and velocities from the tree. As tenter is defined as the time when a satellite galaxy passes the virial radius of the main halo, the initial distance is always equal to the virial radius at tenter.
However, we have to divide between those galaxies that are already within the main halo at zi and those galaxies that enter the halo at a later time. For those galaxies that have entered the halo before zi, we expect that they have already lost some of their angular momentum due to dynamical friction and are therefore closer to the central galaxy than the virial radius. Therefore we scale their initial distance and velocity by where ti is the starting time of the simulation, tenter is the time the satellite entered the main halo and t df is the dynamical friction time (i.e. the time it takes the satellite to merge with the central galaxy). The direction of the position and velocity vectors are unaltered. Before the simulation is started we include all satellite galaxies that have entered the main halo before ti and all satellites that are entering the halo at ti in the initial conditions. We evolve this system with the hydrodynamical code until the time when the first satellite enters the main halo tenter,1 as specified by the merger tree. At this time we interrupt the evolution of the simulation and include the particle realisation of this satellite galaxy in the simulation. The simulation is then resumed and this process is repeated until the final time of the simulation t f . The main advantage of including the satellite galaxies only when they enter the main halo is that the computational cost is reduced. The satellites are only simulated when they affect (or are affected by) the main system. Up to this point we use the SAM to compute the evolution of the satellite galaxy progenitor's properties. In Table 1, we summarise the free parameters of our model and present our fiducial values. We have tested these values mainly on merger trees of MW-like galaxies. For systems of higher or lower halo mass they may have to be adjusted accordingly. We also list the free parameters that are used in the hydrodynamical code and its cooling, star formation and wind models.

SIMULATIONS OF MILKY WAY-LIKE GALAXIES
We employ our model in a series of simulations of MWlike galaxies, i.e. disc galaxy systems with a dark matter halo mass of M h ≈ 10 12 M . The aim of this study is to investigate how a central disc galaxy evolves from redshift z = 1 to z = 0 while experiencing minor mergers. We use dark matter merger trees drawn from the Nbody simulation presented in section 3. As we are interested in systems that experience only minor mergers after z = 1, we select 10 trees with a final halo mass of M h ≈ 10 12 M , that have only mergers with a mass ratio of µ < 0.2. The properties of these trees are presented in Table 2. As a starting redshift we use zi = 1 which corresponds to a cosmic time of t = 6.25 Gyr. We employ the fiducial model parameters presented in Table 1 and simulate all trees up to a redshift z = 0.
The surface density for a typical merger tree (2126 ) is shown in Figures 7 (face-on projection) and 8 (edge-on projection) for the stellar component (upper panels) and the gaseous component (lower panels). The central galaxy evolves in isolation for ∼ 1 Gyr until the first two satellites enter the halo. These satellites pass the central galaxy at a pericentric distance of 55 and 90 kpc, respectively. While the first satellite continues orbiting around the central galaxy and subsequently loses angular momentum and energy, which leads to a decrease of the pericentric distance, the second satellite has enough orbital energy to leave the main halo again. Near the end of the simulation the last three satellites enter the halo and pass the central galaxy.
During the evolution of the simulation, the stellar disc of the central galaxy is heated and thickened. At the same time, a stellar halo forms as a result of two processes. First, due to the close passage of the satellites, some stars in the disc are scattered out of the disc, as energy and angular momentum is transferred. Second, some satellite stars are accreted onto the central galaxy. These stars still retain some of the angular momentum of the satellite such that they do not settle in the disc, but remain in the halo. Additionally, we can identify several stellar streams which originate from the destruction of satellites. However, due to the limited resolution in this pilot study, the number of particles in these streams is not very high, so that a detailed study is not possible. This limitation can be overcome in future studies by increasing the parameter Nres,sat, so that there are more particles in the satellites and thus more particles in the streams. When satellite galaxies enter the halo, a fraction of their cold gas is stripped and eventually can fall towards the centre of the system. Thus some of the stars in the central galaxy form from gas that has been stripped from satellites and accreted onto the central gaseous disc.   Table 2. Properties of the galaxies in the merger trees. For the main system the properties are given at the starting time of the simulation t i , while for the satellites the properties are given at the time when they enter the main halo tenter. All masses are in solar units and all scales are given in kpc.

Evolution of the stellar mass
In this pilot study, we analyse the properties of the central galaxy in each merger tree and these galaxies' evolution during the course of the simulation. Clearly, the stellar mass of the central galaxy is an important quantity. In order to measure the stellar mass in the simulation we identify FOF groups for every snapshot with a linking parameter of b = 0.2. Each FOF group is then split with the Subfind code into a set of disjoint subhaloes that contain all particles bound to a local overdensity. The main halo which contains the central galaxy is the most massive subhalo of the FOF group. The stellar mass of the central galaxy is the total mass of all stellar particles that belong to the main halo. We plot the resulting evolution of the stellar mass for all merger trees in Figure 9. The results of the simulations are given by the red line, while the semi-analytic prediction is given by the blue line. We find that all systems increase their stellar mass and have a final stellar mass that is more than twice the value at z = 1. The stellar masses found in the simulation at z = 0 range from log(m * /M ) = 10.3 to 10.7, while the majority of systems have a stellar mass close to log(m * /M ) = 10.5. These results are in excellent agreement with those obtained by Moster et al. (2010a) using a statistical halo occupation model. For some merger trees the prediction by the SAM agrees very well with the result of the simulation (e.g. the trees 1775, 1975 and 2126 ), while for other trees the semi-analytic prediction for the final stellar mass is much higher than what we find in the merger simulation (e.g. the trees 1968, 1990 and 2048 ). We note that the semi-analytic values for the stellar mass tend to exceed, by a small amount, the average observed stellar mass for a halo of mass 10 12 M . This is because the version of the SAM that we used had been optimized for use with different merger trees, and we did not re-tune the parameters to specifically match this constraint.
In the evolution of the stellar mass, we see that at some times the mass increases quickly in a short time interval. This is not due to an increased SFR due to an interactiontriggered starburst, but is a result of the accretion of the satellite itself. The mass of the accreted satellite is much larger than the newly formed mass as a result of a starburst. Examples can be seen for tree 1975 (t = 11.5 Gyr), tree 2048 (t = 9 Gyr) and tree 2181 (t = 9.6, 11.2 and 12.1 Gyr). Overall however, we find that most stars in the central galaxy have formed from the cold gaseous disc that is accreted from the gaseous halo (in-situ formation), and a only few stars in the central galaxy originate from accreted satellites (ex-situ formation).

Evolution of the cold gas fraction
We further study the evolution of the cold gas component, i.e. all gas particles with a temperature below T = 2 × 10 4 K that reside in the central disc. The resulting gas fractions are plotted in Figure 10. We find that all systems have a higher gas fraction at z = 1, with typical values of 30 to 50 per cent. Towards z = 0 the gas fraction then decreases to lower values of 20 to 40 per cent. These results are in very good agreement with those obtained by Stewart et al. (2009) based on indirect gas fraction estimates from star formation rate densities.
In all systems that have just accreted a subhalo, the gas fraction is elevated. This happens during the first passage of the satellite, in which a considerable amount of its cold gas is unbound by tidal stripping. This gas then falls onto the central galaxy, leading to an increased SFR. We note that unlike in the SAM, this usually happens during the first passage, and not in the final coalescence, which can be much later. The semi-analytic predictions for the gas fractions agree very well with those that are found in the simulations. This is due to the very similar cooling rates in the halo for the SAM and the simulations. The accretion of cooled gas from the halo is much larger than the accretion of stripped cold gas from satellite galaxies. Although the amount of accreted gas from satellites is different in the SAM and the simulations, this effect is small, so that the total gas fractions in the central discs agree quite well.

Evolution of the scale length and height
Finally we study the evolution of the scale parameters of the stellar components, i.e. the scalelength and height of the stellar disc. In order to measure these quantities in our simulations, we first use the decomposition of the particles into disc and spheroidal components. For the disc component we then compute the face-on projected surface density profile and fit the exponential disc scalelength. Similarly, we compute the projected edge-on surface density at a radius of R = 8 kpc, and fit a sech 2 function to it to obtain the disc scaleheight.
We plot the resulting scale parameters in Figure 11 for all merger trees. The disc scalelength R disc is given by the red line and the disc scaleheight z0 is given by the green line. For comparison we included the prediction for the disc scalelength by the SAM which is shown by the blue line. For some systems, the scalelength in the simulation and the SAM agree very well (trees 1775, 1990 and 2092 ), while for other systems the scalelength predicted by the SAM differs from the one found in the simulation (trees 1808, 1877 and 1975 ). Interestingly, while the scalelength in the SAM always increases with time, the scalelength measured in the simulation can both increase and decrease. Overall, we find a broad range of scalelengths for MW-like systems, consistent with observations ( Barden et al. 2005).
The disc scaleheight evolves slowly with time, as long as no satellite merges with the central galaxy. For most systems the ratio between the scaleheight and scalelength is roughly constant throughout the simulation. Final scaleheights range from z0 = 0.6 to 1.2 kpc, consistent with results form observations of MW-like systems (Schwarzkopf & Dettmar 2000;Yoachim & Dalcanton 2006). In some merger trees (trees 1975 and 2048 ), the thin disc is completely destroyed due to mergers and only a thick disc with a final scaleheight of ≈ 4 kpc remains.

DISCUSSION AND OUTLOOK
In order to study the connection between the internal structure of galaxies and their formation histories and large scale environment, a galaxy formation model must have high resolution, include gas physics, and include the cosmological background. The purpose of this paper is to present a novel approach to study the evolution of galaxies by combining semi-analytic models with numerical hydrodynamic merger simulations. Using the predictions of the semianalytic model as the initial conditions for multiple merger simulations, we were able to achieve high resolution at a fraction of the computational cost of standard cosmological hydrodynamic methods.
Hydrodynamical simulations of binary galaxy mergers have been used extensively to study the evolution and morphological transformation of single galaxies due to merger effects, although additional information must be used to place them in a cosmological context. They can be regarded as a method to study the effects of a single encounter in a detailed manner, rather than the complete evolution of a galaxy. For these merger simulations one creates pre-formed galaxies, that are stable in isolation, with parameters drawn from a grid or motivated by observations. These model galaxies are then set on an orbit and evolved with a hydrodynamical code.  Although this technique is useful for studying the evolution of a galaxy during a single encounter, it is not able to predict the typical galaxy properties at a given redshift, as the cosmological background, i.e. the merger history of each galaxy is not taken into account. For this reason, we studied the merger histories of galaxies using merger trees drawn from a large cosmological N -body simulation. We studied whether most mergers are binary or whether they usually involve multiple galaxies, i.e. if two galaxies generally have enough time to merge after their haloes have merged, before the remnant merges with another galaxy. We found that the probability for two mergers happening within a few halo dynamical times is always higher than the probability for many dynamical times to elapse between mergers. In our merger trees, for more than 50 per cent of all merger pairs, the second satellite enters the halo within three dynamical times after the first satellite entered, independent of mass ratio. This indicates that multiple mergers are more common than sequences of isolated binary mergers. As a consequence, it is not sensible to focus on binary mergers, but rather on merger simulations which consider several satellites that enter the parent halo one after each other.
Due to the large number of parameters involved in mergers (orbit, masses, gas fractions, merger ratios, etc.) it is impossible to cover the whole parameter space in merger simulations by drawing the parameters uniformly from a grid. Instead, a simple and elegant path is to use semianalytic galaxy merger trees to generate the initial conditions for galaxy merger simulations. In this way, we automatically select mergers that are expected to be common in the Universe. However, in order to do this, we had to extend the code that creates the particle representations of galaxies used in simulations, as smooth accretion of dark matter and the hot gaseous halo component were not taken into account before. In order to model the cooling and accretion of gas onto the disc we included a slowly rotating hot gaseous halo in the initial conditions generator. We modelled the smooth accretion of dark matter material that is too small to be resolved as a halo in the merger trees, by placing additional dark matter particles around the halo. These particles were placed into small spherical systems with a Gaussian density profile to represent the many sub-resolution systems that are expected to be accreted. The distance to the halo centre was chosen such that they fall into the virial radius of the halo at a specified time extracted from the merger tree. In order to model a system that is stable in isolation, we computed the distribution function for a Gaussian profile and the velocity of each particle with a rejection sampling technique. We tested this smooth accretion model for an isolated MWsized halo and found an excellent agreement with the results from a full cosmological simulation. Furthermore we verified that the small spherical systems are quickly dissolved once they enter the halo.
With the extended model for the initial galaxies, we were able to develop our novel approach that uses semianalytic predictions as initial conditions for a multiple merger simulation. Choosing a starting redshift, we first created particle representations of the central galaxy at this time, and the satellite galaxies at the time when they enter the main halo. We set the mass resolution by requiring that the final number of stellar particles equals a fixed parameter, provided the final stellar mass in the SAM and the simulation are equal. The multiple merger simulation was then performed by evolving the central galaxy until the first satellite enters the main halo, at which point the particle realisation of the satellite was included in the simulation at the virial radius. After this the simulation was resumed and the procedure repeated until the end of the run. We applied our method for ten merger trees drawn from a large dissipationless N-body simulation, for which the main systems are MW-sized at z = 0 and have no major merger after the starting redshift z = 1. For all ten trees we analysed the evolution of the properties of the central galaxy and found good agreement between the stellar mass in our simulations at z = 0 and observational constraints. Similarly, the disc scale lengths in our simulations agree well with the observed range of scale lengths. Overall we also found good agreement with the predictions by the SAM.
There are many differences in the details of the star formation and stellar feedback prescriptions in the hydrodynamic simulations and the SAMs, and the values of the parameters controlling the subgrid physics, and we made no attempt to tune these parameters to match. The fact that there is none-the-less reasonable agreement between the important parameters of the galaxies indicates that the two tools can be used together despite the differences in the details. In the future, it would be straightforward to modify the physical ingredients of the SAM and/or hydro simulations to achieve even better agreement.
Although we demonstrated the basic elements of our new approach here, there are still a number of improvements that can be made and physical processes that can be included. For example, cosmological hydrodynamic simulations have shown that filamentary large scale structures can efficiently feed cold gas into massive halos via "cold streams" (Kereš et al. 2005;Dekel & Birnboim 2006). Because our method assumes accretion from a spherically symmetric halo, it does not include this accretion mode. This process is dominant in massive haloes at a redshift of z > 2. Thus one possible improvement of our method is to model these cold streams.
Another improvement is the inclusion of additional or more realistic feedback processes. In our simulations we have implemented winds with constant mass loading factor and constant wind velocities for all galaxies. It has been shown that better agreement with observations can be obtained with more sophisticated wind models, in which the mass loading factor and wind velocity scales with galaxy properties, e.g. momentum-driven winds (Oppenheimer & Davé 2006). More detailed chemical evolution recipes that follow multiple elements can be implemented to track abundance ratios in various galactic components. Another form of feedback that is not considered in our simulations yet is AGN feedback. The processes that result from gas accretion onto a central BH can provide a large amount of energy, which may be able to reduce cooling from the hot halo after a major merger.
We plan to use the results of simulations carried out with our new approach to study the efficiency of starbursts and morphological transformation (formation of spheroids) in mergers. This will lead to more physical recipes that can then be fed back into the SAM.
In conclusion, we anticipate that the approach developed here will open a path to many possible future applications, including (but not limited to) the study of the formation of stellar halos through the accretion and stripping of satellites, sub-structure and streams in galaxy halos, the build-up of spheroids and thickening of galactic discs through mergers, and the evolution of spheroid sizes and internal velocities through mergers.