Analytic and numerical realisations of a disk galaxy
Recent focus on the importance of cold, unshocked gas accretion in galaxy formation – not explicitly included in semi-analytic studies – motivates the following detailed comparison between two inherently different modelling techniques: direct hydrodynamical simulation and semi-analytic modelling. By analysing the physical assumptions built into the gasoline simulation, formulae for the emergent behaviour are derived which allow immediate and accurate translation of these assumptions to the galform semi-analytic model. The simulated halo merger history is then extracted and evolved using these equivalent equations, predicting a strikingly similar galactic system. This exercise demonstrates that it is the initial conditions and physical assumptions which are responsible for the predicted evolution, not the choice of modelling technique. On this level playing field, a previously published galform model is applied (including additional physics such as chemical enrichment and feedback from active galactic nuclei) which leads to starkly different predictions.
Two different approaches are widely used to test theories of galaxy formation. Both make use of developing computational resources to integrate a set of coupled differential equations forward in time, where each equation applies physical constraints, or empirical laws, to some of the system’s properties.
The semi-analytic method (?) can be characterised by the intent to encapsulate as much of the physical behaviour as possible in the equations themselves, chosen to describe systematic properties (virialisation, conservation of angular momentum, etc.) and thus leaving the minimum to numerical calculation.
Inherently different to this is the approach taken by simulations, where the ethos is to apply equations which are known to govern the particular properties (e.g. the gravitational force between two particles) and demonstrate, by skilfully programmed numerical calculation, that these reproduce macroscopic properties.
In practice, adjustment will always need to be made to account for the fact that each simulated particle represents many, many real particles. In traditional simulations which considered only gravity, this adjustment is purely mathematical. When gravity is the dominant interaction, results are therefore very robust and simulations are mostly accepted as being a faithful realisation of the Universe at large scales, where this is the case. A notable measure of this acceptance is that semi-analytic models have taken to using samples of halos taken directly from simulations (?; ?; ?), rather than samples created statistically (?; ?) using methods based on the analytic arguments of ?.
In general, though, there must come a scale below which the governing equations are no longer describing particle interactions but emergent phenomena. At these smaller scales, the two largely contrasting approaches begin to look more similar. As an immediate example, simulations which encompass an entire galaxy have not attempted to resolve individual stars and must therefore apply an equation which approximates the conversion rate of large bodies of gas into stars, based on locally averaged density. An equivalent and sometimes identical equation appears in semi-analytic models (§3.7).
Both approaches have had some success in accounting for the observed properties of galaxies. Predictive power, on the other hand, relies on surety in the underlying equations, or a clear understanding of how differing underlying assumptions map on to final observable quantities.
To bolster this crucial understanding, it is extremely helpful to study, in detail, the predicted evolution of a single, large galaxy and its satellites, as followed by these two different techniques. If the analytic model begins from the precise halo merger history found in the simulation, and applies the same assumptions about the physics at smaller scales, will it predict the same observable properties and, if not, where and why do the two calculations diverge?
1.1 Previous research
Existing quantitative comparisons between simulations and semi-analytic models have focussed on samples of galaxies. The approach which is adopted here – adapting the semi analytic model to emulate the assumptions made by the simulation—was also used by ? to demonstrate that the cooled mass predicted by the two approaches can be consistent. In particular, consistency was demonstrated over two orders of magnitude in halo mass. Earlier works by ? and ?, using similar techniques, reached broadly similar conclusions.
A broader comparison was made by ? of two predicted galaxy populations: one by an SPH model and the other by version of the GalICS semi-analytic model, modified to emulate the SPH method. Again, consistency is demonstrated between the two techniques once the semi-analytic model is appropriately “stripped down”, but this consistency comes at the expense of agreement with observational constraints, which has been met by the full version of GalICS.
To extend this previous research, in this present work we will compare in detail the formation of a single galaxy in a gasoline numerical simulation and in the galform semi-analytic model. The structure of this paper is as follows:
§2: A review of the broader aspects of the gasoline code and of this particular simulation.
§3: Detailed analysis of the way that gasoline models each of the following key physical processes, together with an explaination of how the particular approach in question can be emulated in the galform model.
§3.3: UV heating
§3.4: Gas distribution
§3.6: Disk formation
§3.7: Star formation
§3.9: Disk stability
Particularly relevant aspects of the galform model are reviewed (or referenced) when appropriate. General aspects of the model can be found in ? and ?.
§4.1: By applying these same underlying physical assumptions to galform, the two modelling techniques are shown to make broadly consistent predictions.
§4.2: Having demonstrated this consistency, galform is then used to find how this halo would evolve under the physical assumptions made in previous studies (which have been shown to successfully match the collective properties of galaxies).
The “MW1” simulation analyzed in this paper was generated using Gasoline, an N–Body+Smoothed Particle Hydrodynamics (SPH) Parallel Treecode (?), within a flat, -dominated WMAP1 cosmology (, =0.7, , , ). A moderate resolution version of this run was originally discussed in detail in ?, and the higher resolution version analised here has subsequently been used for a number of studies (?; ?). High mass and force resolution (see Table 1) was achieved using the volume renormalization technique (?).
This simulated galaxy matches the observed mass–metallicity relationship at varying redshifts (?; ?). This success demonstrates that these simulations overcome the historic “overcooling problem” (?), in which gas in simulations cools rapidly, forming stars too quickly and early, and thus producing a mass–metallicity relation that is too enriched at a given stellar mass. ? showed that reproducing the mass–metallicity relationship was a result of the adopted feedback prescription. This feedback mechanism, combined with high numerical resolution111With total particles within the virial radius at z=0., avoids spurious loss of angular momentum. As a consequence, simulated galaxies also match the observed Tully-Fisher relationship (?; ?; ?).
|Dark Matter particles:|
|Force resolution:||0.3 kpc|
Gasoline’s star formation and supernovae feedback scheme is described in detail in ?, and parameters from ? are adopted for MW1. The star formation recipe reproduces the Kennicutt-Schmidt law, and each star particle is born with a Kroupa initial mass function (?). Combined with a prescription for a cosmic UV background (?) that mimics reionization and unbinds non-collapsed baryons from simulated halos with total masses below about , this feedback scheme drastically reduces the number of luminous satellites containing a significant mass of stars, avoiding the well known “substructure problem” (?; ?; ?; ?).
As described in ?, supernova energy is deposited into the ISM mimicking the blast-wave phase of a supernova (following the Sedov-Taylor solution). Cooling is turned off for gas particles within the supernova blastwave radius for the duration that the remnant is expected to adiabatically expand (§3.8). Supernova feedback regulates star formation efficiency as a function of halo mass, resulting in the mass–metallicity relationship described above. This regulation of star formation also leads to realistic trends in gas fractions, with the lowest mass galaxies being the most gas rich (?), and reproducing the observed incidence rate and column densities of Dampled Lyman systems at 3 (?).
MW1 and its most massive progenitor halo were identified at each simulation output step using AHF222AMIGA’s Halo Finder, available for download at http://popia.ft.uam.es/AMIGA (?; ?). AHF determines the virial radius, , for each halo at each output time step using the overdensity criterion for a flat universe following ?. The full gas accretion history of MW1 is described in detail in ?, and used in this study for comparison to Galform.
3 Translation from simulation to semi-analytics
3.1 The merger tree
The common point between the two different models of galaxy formation is the halo merger tree. This is generated by finding bound structures in the “MW1” simulation of ? using AHF. A virial mass is assigned to each structure and its fate at each timestep is categorised, being either:
Inside the virial radius of a larger system; or,
Two criteria are used to classify a galaxy as “merged”:
If the DM mass decreases by more than 50% in a single timestep; or,
If the number of DM particles falls below 64.
These are both significantly different from the definition adopted by galform (see ?), so a detailed comparison of merger times, along the lines of (?), is deferred for another study.
When a structure changes from category 1 to category 2, the galform code updates its status to a “satellite” and predicts its subsequent evolution based on its approach velocity, v and position r at this juncture. These two properties are combined into the dimensionless parameter (§32), which characterises the orbit. The distribution of values for in this simulation was found to be consistent with the most recent assumption applied in galform (Fig. 1).
This orbital information is then used to estimate a merger time based on the standard Chandrasekhar formula (37) as defined by ?. Details of this part of the model can be found in Appendix A and in ? (?; §4.3).
The resulting dark and baryonic components are shown in Fig. 2. Being generated purely from the derived merger tree, the two realisations might be expected to be identical. One reason this is not quite the case is because the total virial mass, , is measured directly, but the relative composition of baryons and dark matter may vary a little from halo to halo, and from one realisation to the other.
Previously, galform required that the halo mass was monotonically increasing as a function of time along each branch of the merger tree (i.e. each halo should be more massive than the sum of the masses of its progenitor halos). While this condition is always fulfilled for merger trees generated using the extended Press-Schechter theory, it is not necessarily fulfilled for merger trees extracted from N-body simulations.
For this reason, previous use of N-body merger trees in galform has required some adjustment of halo masses333The standard approach which was taken begins at the earliest times in the tree, checks whether each halo is more massive than the sum of the masses of its progenitor halos and, if it is not, sets its mass to be very slightly larger than that sum. This process was repeated for all halos in the tree. to enforce this monotonic behaviour (?). Unfortunately, this adjustment is somewhat arbitrary and can lead to significant changes in halo mass when adjustments are propagated through the merger tree.
Monotonicity in halo masses along branches of the merger tree is not enforced in this work. Instead, they are allowed to follow whatever evolution the N-body simulation provides. This requires some modification of the way in which galform assigns baryonic masses to halos. Previously, each halo was assigned an initial mass of hot gas equal to:
where was the total mass of the halo and “progenitors” refers to halos that exist on a previous timestep but have, at some point, merged into the current “parent” (or its most massive progenitor). If halo masses do not monotonically increase, this could lead to negative gas masses. To ensure against this, the initial mass of gas is assigned to be the greater of (1) and zero.
The results are identical to the original implementation in cases where the halo masses are monotonically increasing, but this can result in the total baryonic mass not being precisely equal to the fraction of the total mass of the final halo. Fortunately, this difference is found to be small and, in any case, the simulated halos do not have precisely this ratio either.
The temperature history of accreted gas was investegated by ? in a simulation containing 1120 galaxies. The distribution of gas particles in terms of their maximum past temperature is bimodal, suggesting that classification into hot and cold modes of accretion is appropriate. The contribution from each of these modes is then found to depend strongly on halo mass and on environment. Notably, the halo mass () which separates low-mass galaxies dominated by cold accretion and higher mass galaxies dominated by the hot mode, is close to that found by analytic arguments (?).
The decomposition of the accreting gas into different modes is a principle feature of the analysis by ?. The three modes of accretion are defined as follows (along with their analogues in the semi-analytic realisation):
Clumpy acccretion is defined to include “any particles that have, at any previous output step, belonged to a galaxy halo other than the main halo we are considering.” (?)
The information needed to evaulate this component in the semi-analytic model is readily available from the merger tree, being the virial masses, , of the existing, progenitor halos:
Smooth accretion: Shocked component
The following criteria are required for a gas particle in the simulation to be categorised as shocked:
where is the virial temperature defined in gasoline
and is the fractional increase in density since the previous timestep. The latter is an entropy criterion and the factor of 4 is motivated by the Rankine-Huginot shock jump conditions for a singular isothermal sphere. Though low resolution can artificially broaden shocks, ? verify that this run is of sufficiently high resolution to resolve shock discontinuities.
The virial temperature in galform is defined differently, being a factor of higher, , so (3a) corresponds to .
Finding an appropriate analogue to this condition within the framework of galform is difficult. However, one relevant calulation that is made compares the cooling timescale for gas at temperature :
and the free fall time from a radius, , enclosing mean density :
The cooling function is found for given temperature and metallicity, , by reference to the published tables of ?. is the gas density, the mean density of all the mass enclosed by radius , and is the mean particle mass.
So, in the context of galform accreting gas is designated as shocked if:
this being a reasonable, simple equivalent to (3a).
Smooth accretion: Unshocked component
Gas in the simulation which fails to meet either of the criteria in (3) is categorised as unshocked. Similarly, accreting gas in galform is designated as unshocked if .
The relative contributions of these three components are shown in Fig. 3 for both realisations. Being determined directly from the merger tree, it is unsurprising that the clumpy components are quite consistent.
That the shocked component is overestimated by the semi-analytic model is entirely as anticipated by ?, though the eventual impact this has on the predicted properties of the galaxy may not be quite as severe, as will be seen in the forthcoming sections.
3.2.1 Alternative criterion
? investigate the conditions for the existence of virial shocks in galactic halos, culminating in the following expression for the existence of a shock:
As can be seen from Fig.3, the cooling time does not exceed the free-fall time by a factor of 5 until about 8 Gyr. The criterion (9) would therefore imply that all smoothly accreted gas will be unshocked before that point, corresponding to an unshocked:shocked ratio of 2 or 3, as opposed to in gasoline or according to (7).
So, for this particular mass aggregation history, one simple estimate (7), overestimates the shocked component in the simulation and another, more thorough estimate (8), gives an underestimate. This suggests that the error in assuming spherical symmetry may unfortunately exceed any accuracy that can be gained by analysis of shock conditions. Calibrating some of these analytic models to simulations may help quantitatively account for the complex geometry.
3.3 UV heating
A UV background is incorporated by ? which is based on an updated model by ?. This background “turns on at ” and heats the gas in halos so that only those with “ 100 or more dark matter particles are of sufficient mass to retain bound gas particles”. For the resolution in this case, 100 DM particles corresponds to .
The galform model is build with two corresponding parameters: A redshift, , below which the UV background heating is effective (set to for consistency with the above) and a halo velocity, below which all cooling is suppressed. For a virial overdensity, , the mass cut-off, , mentioned above translates to a cut-off halo velocity of:
These values are duly applied, but one key difference between the modelling of this process remains. The UV background in the simulation can slow cooling in larger halos without necessarily completely suppressing it, while in galform suppression is all or nothing (below and above respectively).
3.4 Hot gas distribution
The distribution of hot gas that is assumed in the galform model was recently updated, in the light of work by ?, to the following density profile for hot gas of mass :
This form was assumed in recent work on disk formation and structure (?), and will also be used here.
The reason for this choice is clear from Fig. 4, which demonstrates that it is a satisfactory description of the distribution of diffuse gas in the outer parts of this simulated system. In this Figure, the analytic form for the hot halo gas (11) and cold disk gas (21) have been fitted to the simulation results, yielding the respective scalelengths:
Halo core radius, .
Disk scalelength, .
(with the initial constraint on each line that it integrates to the correct total mass for that component.) These scalelengths, found from fitting the simulation results, can be scrutinised against those deduced from the usual analytic assumptions.
Firstly, the halo core radius is typically assumed to be a constant fraction of the virial radius, which is not quite the case in this particular simulation (Fig. 4). However, in the absence of any obvious pattern in the evolution of , the usual simple scaling is applied, with the choice:
being the best approximation. Fig. 5 compares this to values found from fitting to the true gas distributions in Fig. 4, the result looking favourable. This straightforward scaling with the virial radius (12) will therefore be used in the matched galform realisation. The cold gas scalelengths, also shown in Fig.5, are not inputs to galform, but are calculated within the model as described in §3.6.
To understand the relevance of the hot gas distribution on the system’s evolution as a whole, it is worth reviewing the precise treatment of gas cooling in galform.
The cooling model described by ? determines the mass of gas able to cool in any timestep by following the propagation of the cooling radius in a notional hot gas density profile444We refer to this as a “notional” profile since it is taken to represent the profile before any cooling can occur. Once some cooling occurs presumably the actual profile adjusts in some way to respond to this and so will no longer look like the notional profile, even outside of the cooling radius.. This profile is fixed when a halo is flagged as “forming” and is only updated when the halo undergoes another formation event. The mass of gas able to cool in any given timestep is equal to the mass of gas in this notional profile between the cooling radius at the present step and that at the previous step. The cooling time is assumed to be the time since the formation event of the halo. Any gas which is reheated into or accreted by the halo is ignored until the next formation event, at which point it is added to the hot gas profile of the newly formed halo. The notional profile is constructed using the properties (e.g. scale radius, virial temperature etc.) of the halo at the formation event and retains a fixed metallicity throughout, corresponding to the metallicity of the hot gas in the halo at the formation event.
This work makes use of a new cooling model, which will be described in full detail in ?. Rather than arbitrary “formation” events, this model uses a continuously updating estimate of cooling time and halo properties. The properties described in §3.4 (density normalization, core radius) are reset at each timestep. The previous infall radius555defined as the lesser of the cooling and freefall radii. (i.e. the radius within which gas was allowed to infall and accrete onto the galaxy) is computed by finding the radius which encloses the mass previously removed from the hot component in the current notional profile.
The new model must supply an alternative estimate of the time available for cooling in the halo, , from which the cooling radius can be computed in the usual way (i.e. by finding the radius in the notional profile at which ). This is done by considering the energy radiation rate per particle, and the thermal energy per particle, , which are estimated by making standard assumptions:
In terms of these quantities, the cooling time (5) is simply:
At any time, the model needs to identify some radius in the hot halo, , where the gas has had just enough time to radiate all its thermal energy. This as yet undetermined radius is defined to satisfy the condition that the cooling time for gas particles at this point is equal to the time available for them to cool, which is estimated in terms of the mean energy radiation rate per particle666This is proportional to the mean number density : (15) The factor does depend on halo structure (i.e. on ) but, if the distribution of hot gas relative to the virial radius is the same for each halo, as implied in this case by (12), then this factor cancels out in (17)., :
where is the mean number density inside the virial radius.
This approach must account for the hierarchical assembly of a halo, so the denominator in the right hand side strictly includes a sum over all progenitors, , of this system which exist at time . Incorporating this sum into (17), and rearranging, then yields a final expression to be solved for the cut-off density, :
This equation is notably independent of halo structure, a factor which only comes in when the limiting density is used to compute the total cooled mass. The scalelength, , appears in the normalisation of the hot gas profile (11), and hence effects the mass enclosed by the cut-off density, .
This can be appreciated graphically. Solutions to (18) are illustrated in Fig. 6 for three choices of and three examples of cut-off density. For the higher cut-off density, a smaller (more centrally concentrated) will lead to a greater cooled mass. For the case of a low cut-off density, a centrally concentrated profile will actually lead to a lesser cooled mass, contrary, perhaps, t
3.6 Disk formation
In galform, galactic disks are assumed to form from the cooled gas, conserving the net angular momentum that it possessed when it was part of the halo. The assumed distribution for the specific angular momentum, , of the hot halo gas is taken from ?:
The value of was found to lie between 0.5 and 1.5, but the mean value is adopted here. This is normalised to enclose the correct total mass and the mean specific angular momentum of the hot gas, , which is found directly from the spin parameter,
Details of the calculation of the halo energy, , can be found in ? (?; §3.2.1). The disk scalelength, , is recalculated after each timestep such that angular momentum is conserved in the system. The calculation is based on the following surface density distribution (applied to each disk component, ):
and a circular speed, , which is constant with radius. The subsequent evolution of the disk is then followed by applying physical assumptions which are as close as possible to those made in the simulation, as is explained in the following sections.
As with the sub-halo trajectories, the parameter can ordinarily be assigned at random, from a specified distribution. In this case, its value is measured for each of the simulated halos and passed on to galform, along with the rest of the merger tree information (§3.1).
The goal here is not to prove that parameters like , and must necesarily be set by simulations in order to validate the results of any semi-analytic model. Rather, it is an example of the correspondence between these initial conditions and the properties of the galaxy which they produce. This should reassure theorists that, when such models are applied to large samples of galaxies, their individual properties will be representative as long as these governing, random, but physical parameters are drawn from the correct natural distributions.
3.7 Star formation
At face value, gasoline and galform would appear to be applying the same assumptions regarding the conversion of cold gas to stars. These respective assumptions can be found from equation (4) of ? and equations (4.4) and (4.14) of ?777An additional factor involving the circular speed appears in (4.11) of ?. This is omitted here (equivalent to setting ).:
where is the total instantanoues star formation rate in the whole disk, and is the local value per unit volume. The timescales which appear in these expressions have similar physical motivation, but are not identical. One, , is the time for gravitational collapse of a spherically symmetric region of local density . The other, , is the angular period of the whole system at the half-mass radius, :
So, despite the apparent similarity of the two expressions in (22), they can lead to quite different instantaneous total star formation rates, . To illustrate this, one can consider the surface density distribution (21) which galform assumes, together with some decrease in density away from the plane of the disk:
Applying the local star formation assumption of (22a) to this distribution gives:
Comparison between this and the usual prescription in galform exposes the big physical difference between the two star formation models; (22a) considers the gravitational collapse of the gas due to its own gravity while (22b) includes the gravity of all components888The dark matter does influence the creation of stars due to its effect on the overall stability, but this is a separate part of the analytic treatment (§3.9)., including the stars and dark matter ( rather than just ).
The derived analytic equivalent (25) is shown in Fig.7 (using the fitted values of from §3.4 and the illustrative approximation999This choice is not too critical since appears only to the one half power in eqn. 25). of a constant scaleheight, pc) alongside the full simulation results. Discrepancies are larger at early times when the assumed cold gas gas density distribution is tenuous, but the match improves as the disk settles and is excellent at late times, showing that (25) is as consistent with the assumptions in Gasoline as is possible within this particular analytic framework.
3.7.1 Burst star formation
During merger events, and if the disk is deemed to be unstable (§3.9), cold disk gas in the galform model is assumed to funnel to the centre of the new system and be converted to stars within a timestep, which in this case is about 0.2Gyr (? ?; §4.3). The incorporation of this non-equilibrium effect should broadly correct the total star formation rate on occasions when the assumed density profile (21) is not such a good description of the true gas distribution. This is relevant here in the early stages of the simulation (as apparent when comparing the points and solid line in Fig. 7).
3.7.2 Schmidt-Kennicutt Law
Equation (22b) can be recognised as a version of the established empirical relation between star formation rate and gas supply (?), for which the constant of proportionality has the value . The corresponding star formation efficiency for this simulated system can be calculated from its approximate half-mass radius and circular velocity (Fig. 10) and this is shown in Fig. 7 as a guide.
3.8.1 Heating of disk gas
The MW1 simulation treats each type II supernova event as a blast wave which expands to a maximum radius, , and the enclosed volume is assumed to be unable to cool for as long as the shell can be expected to survive, . If the star formation rate varies little over this timescale, and the shell subsequently cools rapidly, this will correspond to a heated mass,
where is the average mass of stars formed per supernova. Now, (26) conceals the fact that both and (and hence ) are, rightly, functions of the local gas density, and pressure, . The precise dependence assumed is found in equations 9 and 11 of ?. Combining this into (26) gives:
Referring to the phase diagram in figure 4 of ?, gas in star forming regions appears to be approximately isothermal at K (which is expected due to the steep gradient in the cooling curve at this temperature). This leads to the suggestion that the environmental dependence in (27) will be minimal, the pressure and density terms almost cancelling for constant temperature. Though tenuous, this provides at least a rough explaination for the extremely simple dependence which does indeed emerge from the full calculation, as shown in Fig. 8:
The precise definition of used in the figure is the mass of disk gas which fails to satisfy the temperature criterion for star formation (K). Thus, the integrated effect of feedback will be to slightly reduce the fraction of the total disk gas which is available for star formation101010For a massive galaxy such as this, the star-forming timescale is typically years, so the correction to allow for re-heating actually makes very little difference. Quick examination of the inset to Fig. 8 confirms this; the fraction of the total gas which gets heated above the star-forming limit is tiny..
where is the particular timescale derived from equation (25). This revision is duly applied in galform, being a fair reflection of the physical assumptions in the simulation.
To discover the extent of outflow from the galaxy in the simulation, all particles were identified which have moved from the inside to the outside of a radius of 30 comoving kpc. The total mass of such material turned out to be even less than the heated disk mass (see Fig. 8) and can safely be ignored for the purposes of this study. The outflow of gas as a result of supernovae in the central galaxy is therefor set to zero in this version of galform, for consistency with the Gasoline simulation.
It is important to note that this feedback model is more effective at lower halo masses, due to shallower potential wells (?). This dependence is also accounted for in the usual galform feedback model, where the rate of mass outflow is assumed to be at least equal to the star formation rate, and thousands of times higher for smaller systems (?). The consequence of this difference is investigated in §4.3.
3.9 Disk stability
The analytic condition which is used to assess the stability of disks in the galform model follows ?. It deems disks to be unstable to bar formation if:
When this inequality is met, all disk material is transferred to a central bulge.
In the usual galform model, some set fraction of disk mass is also transferred to a central black hole (?). This was used to arrive at the value in (30) which was found to best match the Magorrian relation between bulge mass and black hole mass, , observed by ?. In this study, the black hole accretion fraction is set to zero for the purposes of consitency with the simulation, where this aspect of evolution is not considered.
Being physically well motivated, the condition (30) has been applied throughout this study. Conveniently, the precise formalism used need not cause too much concern here; the change that its application produces in the component masses (in Fig. 11, for example) is barely enough to be noticed. This can be understood from Fig. 9, which shows that only a very limited adjustment to the disk mass is required to prevent the limit in (30) from being reached. The mass aggregation of this particular system, and the fact that the specific angular momentum of the gas is relatively high111111 The criterion (30) can also be written: (31) where and are the masses of the dark matter and spheroid components. The disk’s eventual half mass radius will enclose more dark matter if the original spin parameter is higher. ( for the latter half of the halo’s history) leads to the stabilisation of the system, even if the redistribution of mass is not enforced.
4.1 Matched Model
The size and rotation speed are shown for both realisations in Fig. 10. The predicted galform value is based on conserving the initial angular momentum (19) of in-falling gas; an assumption which is particularly effective in this case where the distributions for the initial and final mass were appropriate (Fig.4). This correspondence is encouraging for the approach of calculating disk sizes from the spin parameter, , of the halo from which the gas cooled.
4.1.2 Cooled Mass
The total cooled mass (stars and cold gas), calculated using the methods of §3.2, can be seen from Fig. 11. Agreement with the simulation is excellent, given the inherent difference in the methods of calculation, and bolsters confidence in the rather intricate calculation of the energy radiated by the system throughout its complex merger history (13-18).
The division of the cooled mass into stars and gas are also shown in Fig. 11. The evident agreement can be understood by reference to Fig’s 4, 7 and 8, which demonstrate that the considerable complexities of the simulation reduce to relatively simple relationships when integrated over the entire disk.
There is a minor discrepancy at early times (Gyr) in the cooled mass generated by the two modelling techniques. This is unlikely to be due to differences in cooling from the hot halo: given the similar density distributions (Fig. 4) and cooling function, the two models should yield similar amounts of cooled gas, all else being equal.
The discrepancy is also not due to a difference in the time of development of a hot gas halo. The top panel of Fig. 4 shows that accreted gas in the galform model transitions from unshocked, cold gas to shocked, hot gas at an age of about 2 Gyr, while the bottom panel shows that this time is also when the gasoline model begins to develop a hot gas halo.
This demonstrates that the adopted analytic separation of unshocked and shocked gas (free fall limited regime versus cooling limited regime) is an excellent analog to the transition between cold gas accretion and shocked gas accretion seen in simulations as a galaxy crosses the threshold mass capable of developing a hot halo. It is important to note that this transition has always been included in analytic models of galaxy formation, and thus cold gas accretion at low galaxy mass does not alter the standard picture of galaxy formation.
The discrepancy in cooled mass arises after the development of this hot halo. The simulated galaxy continues to accrete cold gas via filaments that penetrate within the hot gas halo (see Fig. 3 here, figure 5 of ? and also ?, ? and ?), while cold gas accretion ends in the analytic model.
4.2 Adapting gasoline to emulate semi-analytics
? demonstrated that the inclusion of cold gas accretion along filaments leads to an earlier phase of star formation than predicted if all gas was initially shock heated, due to the shorter cooling times onto the central galaxy of the cold gas. Fig. 12 shows that the simulated galaxy does indeed form stars earlier than the analytic model, though the overall discrepancy is never more than a factor of two.
In an attempt to adapt the simulations to include similar physics as the analytic model, Fig. 12 shows the effect on the stellar growth in the simulations if all of the cold gas had instead been shocked. The hot gas has cooling times several Gyr longer than the cold gas. This is quantified, and a delay has been added to the formation time of the stars spawned from the cold accreted gas (?). As seen in Fig. 12, this delay leads to a slower build of up stellar mass in the simulations. Comparison with Fig.11 shows this to be in even better agreement with the stellar growth of the analytic model, which neglects this cold gas accretion along filaments.
4.3 Comparison with standard galform assumptions
In order to represent the same physical assumptions as this particular simulation, the semi-analytic code had to be significantly modified from the version which has been used most recently to study the collective properties of galaxy samples (?; ?). There are three very significant differences between this version of galform and the simulation studied here, all of which are manifest in the comparison of the history of the simulation and the ? model (Fig. 13).
4.3.1 Inclusion of AGN feedback
When there are instabilities (§3.9) or merger events in this model which trigger disk collapse, 5% of the available disk gas is assumed to accrete onto a central black hole. If the following criteria are satisfied, it is assumed that no hot halo gas will be able to cool onto the disk.
The resulting Eddington luminosity of the black hole is at least 4% of the cooling luminosity of the halo.
This effect is visible in Fig. 13 as a period of rapidly decreasing cold gas mass; during such phases it is no longer being replenished by gas cooling in from the halo.
4.3.2 Stronger supernova feedback
This is assumed to be extremely strong in the ? model, the justification being that a more modest conversion of supernova energy to gas outflow would allow the formation of too many low-mass galaxies. (Indeed, it appears that the predicted number may still be too high (?) even with this generous allocation of feedback energy.)
In conjunction with the inclusion of AGN heating, this assumption leads to a system with a much greater hot gas component: As soon as the supply of gas to the disk can no longer be replenished by cooling, all the gas in the system quickly ends up in the hot component. This effect is additionally enforced by the low star formation efficiency (see below), which prevents disk mass being locked up into stars.
4.3.3 Inefficient star formation
The star formation assumptions in Gasoline and the ? model have already been contrasted in §3.7. The assumed efficiency in the ? model is to be so low () that, even after allowing for the structure and gas fraction, the final conversion rate is about a factor of 10 lower than in the simulation. ( 0 to 0.2/Gyr as opposed to the rate of 0.5 to 3/Gyr seen in Fig. 7).
This stands out in Fig 13. The large mass of cold gas which is supported by the disk in the semi-analytic realisation (except when cooling is shut off by AGN) cannot be sustained in the simulation due to much more effective conversion to stars.
The merger history of a simulated, Milky Way-like galaxy has been used to recreate the same system using the galform semi-analytic model.
The incoming trajectories of satellite halos were investigated and found to be consistent with the usual, assumed distributions. Information on the different modes of accretion onto the main halo is shown to be inherently contained in the merger tree (and would be similarly contained in a statistically generated merger tree which may be used in the study of larger volumes).
The subsequent fate of accreted gas is not modelled equivalently by the two methods, but the analytic framework does account for both shocked and unshocked gas. Furthermore, the estimated formation time of the hot halo is in excellent agreement with the simulation.
The single-parameter analytic distributions which are assumed in galform are found to be a satisfactory description of the radial distribution of gas and stars inside the simulated halo. It was also shown that the choice of hot halo scalelength is of minor importance in the calculation of cooled gas mass.
However, the enforcement of spherical symmetry in these analytic distributions means that any subsequent cold accretion along filaments (persisting in spite of shocks elsewhere) is neglected. Thus, the total quantity of unshocked gas onto the central galaxy is underestimated, as expected by ?. Some gas is consequently delayed in its arrival at the central galaxy, though it is found to proceed at later times.
The disk structure is predicted analytically by conserving the angular momentum of the cooled gas. Resulting scalelengths are within a factor of two of the simulation results and the circular velocities even closer. This is quite satisfactory for the case of a single halo, but more work is needed to establish whether the two techniques will consistently agree to this level for an arbitrary set of initial conditions (and whether any disagreement is biased).
With regard to the subsequent evolution of the disk, the galform code was altered to apply the same fundamental assumptions as the simulation by considering their integrated effect over the assumed disk structure. Subsequent agreement was excellent given this immense reduction in complexity from many thousands of SPH particles to a few one-dimensional equations.
So, by assuming the same physics and using the same initial conditions, the two techniques predict a final system which is recognisably the same, despite fundamental differences in the way that these assumptions are applied and the evolution followed. This suggests that equations which attempt to characterise the emergent behaviour of the system may, with sufficient understanding, become as reliable as those which emulate the properties of particles themselves.
With such potential consistency thus established, the same system was then evolved under physical assumptions which have been more frequently adopted within galform (?). The resulting system (at any time) is not recognisably the same as that predicted by the simulation: having less than a quarter of the stellar mass, barely half the rotational velocity and with negligible star formation at z=0. Though not all aspects of a system need to be correctly represented to learn from the modelling exercise, this divergence from the same initial conditions proves that at least one of these published models poorly represents the true nature of galaxy formation.
It is hoped that this detailed comparison between two very different theoretical approaches will promote the awareness of their respective limitations and help lead them both towards a truer description of real galaxies.
The authors would primarily like to thank the KITP, Santa Barbara for their hospitality, supported in part by the National Science Foundation under Grant No. PHY05-51164, without which this interdisciplinary study would not have taken place.
AMB acknowledges support from the Sherman Fairchild Foundation. AJB acknowledges support from the Gordon & Betty Moore Foundation. FG acknowledges support from the HST GO-1125, NSF AST-0607819 and NASA ATP NNX08AG84G grants.
Thanks go to Richard Bower, Shaun Cole and Carlos Frenk for their very helpful comments. Also, they and their collegues Carlton Baugh, John Helly, Cedric Lacey and Rowena Malbon allowed us to use the galform semi-analytic model of galaxy formation (www.galform.org) in this work, for which we are extremely grateful.
Agertz, O., Teyssier, R., & Moore, B. (2009), MNRAS, 397, 64
Bower, R. G., Coles, P., Frenk, C. S., & White, S. D. M. 1993, ApJ, 405, 403
Bower, R. G., Benson, A. J., Malbon, R., Helly, J. C., Frenk, C. S., Baugh, C. M., Cole, S. & Lacey, C. G. 2006, MNRAS, 370, 645
Benson, A. J., Pearce, F. R., Frenk, C. S., Baugh, C. M., & Jenkins, A. 2001, MNRAS, 320, 261
Benson, A. J. 2005, MNRAS, 358, 551
Benson, A. J. et al. (in prep.)
Birnboim, Y., & Dekel, A. 2003, MNRAS, 345, 349
Blumenthal G. R., Faber S. M., Flores R. & Primack J. R., 1986, ApJ, 301, 27
Brooks, A. M., Governato, F., Booth, C. M., Willman, B., Gardner, J. P., Wadsley, J., Stinson, G., & Quinn, T. 2007, ApJ, 655, L17
Brooks, A. M., Governato, F., Quinn, T., Brook, C. B., & Wadsley, J. 2009, APJ, 694, 396
Cattaneo, A., et al. 2007, MNRAS, 377, 63
Cole S., Lacey C., Baugh C. & Frenk C. 2000, MNRAS, 319, 168
Croton, D. J., et al. 2006, MNRAS, 365, 11
De Lucia, G., & Blaizot, J. 2007, MNRAS, 375, 2
Dekel, A., et al. 2009, Nature, 457, 451
Efstathiou, G., Lake, G., & Negroponte, J. 1982, MNRAS, 199, 1069
Gill, P. D. S., Knebe, A & Gibson, B. K. 2004, MNRAS, 351, 399
Governato, F., Willman, B., Mayer, L., Brooks, A., Stinson, G., Valenzuela, O., Wadsley, J., & Quinn, T. 2007, MNRAS, 374, 1479
Governato, F., Mayer, L., & Brook, C. 2008, Astronomical Society of the Pacific Conference Series, 396, 453
Governato, F., Brook, C. B., Brooks, A. M., Mayer, L., Willman, B., Jonsson, P., Stilp, A. M., Pope, L., Christensen, C., Wadsley, J., Quinn, T. 2008, arXiv:0812.0379
Gill, S. P. D., Knebe, A., & Gibson, B. K. 2004, MNRAS, 351, 399
Gross, M. A. K. 1997, Ph.D. Thesis,
Haardt, F., & Madau, P. 1996, ApJ, 461, 20
Haardt, F., & Madau, P. 2001, Clusters of Galaxies and the High Redshift Universe Observed in X-rays,
Hatton, S., Devriendt, J. E. G., Ninin, S., Bouchet, F. R., Guiderdoni, B., & Vibert, D. 2003, MNRAS, 343, 75
Häring, N., & Rix, H.-W. 2004, ApJ, 604, L89
Helly, J. C., Cole, S., Frenk, C. S., Baugh, C. M., Benson, A., & Lacey, C. 2003, MNRAS, 338, 903
Jiang, C. Y., Jing, Y. P., Faltenbacher, A., Lin, W. P., & Li, C. 2008, ApJ, 675, 1095
Katz, N., & White, S. D. M. 1993, ApJ, 412, 455
Kauffmann, G., White, S. D. M., & Guiderdoni, B. 1993, MNRAS, 264, 201
Kauffmann, G., Colberg, J. M., Diaferio, A., & White, S. D. M. 1999, MNRAS, 303, 188
Kereš, D., Katz, N., Weinberg, D. H., & Davé, R. 2005, MNRAS, 363, 2
Knebe, A. G., Green, A. & Binney J., 2001 MNRAS, 325, 845
Knollmann, S. R., & Knebe, A. 2009, ApJS, 182, 608
Kroupa, P., Tout, C. A., & Gilmore, G. 1993, MNRAS, 262, 545
Kennicutt R. C., 1998, Annual Review of Astronomy & Astrophysics, 36, 189
Maiolino, R., et al. 2008, A& A, 488, 463
Ocvirk, P., Pichon, C., & Teyssier, R. 2008, MNRAS, 390, 1326
Pontzen, A., et al. 2008, MNRAS, 390, 1349
Press, W. H., & Schechter, P. 1974, ApJ, 187, 425
Quinn, T., Katz, N., & Efstathiou, G. 1996, MNRAS, 278, L49
Read, J. I., Mayer, L., Brooks, A. M., Governato, F., & Lake, G. 2009, arXiv:0902.0009
Sutherland R. S. & Dopita M. A., 1993, ApJ, 88, 253
Lacey C. & Cole S. 1993, MNRAS, 262, 627
Maiolino, R., Nagao, T., Grazian, A., Cocchia, F., Marconi, A., Mannucci, F., Cimatti, A., Pipino, A., and 11 co-authors
Mayer, L., Governato, F., & Kaufmann, T. 2008, Advanced Science Letters, 1, 7
Moore, B., Ghigna, S., Governato, F., Lake, G., Quinn, T., Stadel, J., & Tozzi, P. 1999, ApJL, 524, L19
Sharma S. & Steinmetz M., 2005, ApJ, 628, 21
Springel V. 2005, MNRAS, 364, 1105
Stringer M. J. & Benson, A. J. 2007, MNRAS, 382, 641
Stringer, M. J., Benson, A. J., Bundy, K., Ellis, R. S., & Quetin, E. L. 2009, MNRAS, 393, 1127
Stinson, G., Seth, A., Katz, N., Wadsley, J., Governato, F., & Quinn, T. 2006, MNRAS, 373, 1074
Tormen, G., 1997, MNRAS, 290, 411
Wadsley, J. W., Stadel, J., & Quinn, T. 2004, New Astronomy, 9, 137
White, S. D. M., & Frenk, C. S. 1991, ApJ, 379, 52
Yoshida, N., Stoehr, F., Springel, V., & White, S. D. M. 2002, MNRAS, 335, 762
Zolotov, A., Willman, B., Brooks, A. M., Governato, F., Brook, C. B., Hogg, D. W., Quinn, T., & Stinson, G. 2009, arXiv:0904.3333
Appendix A Satellites
The orbits of sub halos around their host can be characterised (?) by the dimensionless parameter
where and are the radius and specific angular momentum of a circular orbit with energy . For a potential with a flat rotation curve; where the circular velocity at all radii is ,
Hence, for a subhalo at position, with velocity , the radius of a circular orbit with the same energy is:
The orbital parameters are found by simply applying equations (32) to (34) to the position and velocity vectors of each satellite halo (after moving to a frame where these vectors are both 0 for the host halo):
The galform model assigns each satellite a value for this parameter, , drawn at random from the log-normal distribution that it was found to follow in simulations by ?:
with . The sample of orbit parameters produced by this random assignment is shown in Fig. 14, together with the actual values from the ? simulation.
The fraction of the halo free-fall time that it takes for the satellite to merge through dynamical friction is assumed to be proportional to :
This is based on the standard Chandrasekhar formula for the dynamical friction and appears originally in ?. is the mean density at the virial radius (determined by the cosmology, not the specific halo’s properties).
The distribution of angles at which sub structures enter their host halo has been studied by ? for the VLS and VIRGO simulations and and found to have a repeatable distribution. This distribution can be applied in the galform model to generate an alternative set of orbital parameters, , which are shown in Fig. 14 alongside the standard assumption of the log-normal distribution, and the results from gasoline.
Merger times predicted by (37) have been compared by ? with the results of GADGET2 simulations (?). The agreement found was of the order of a factor of two, which is deemed acceptable for the continued use of this formula in galform. However, since (37) was found to consistently underestimate the simulated merger time, the improved fitting formula proposed by? may well be adopted in future.
Unfortunately, due to the very different definition of merging adopted by the halo finder used here (§3.1), a meaningful comparison of the respective times in the two realisation has not been possible.
Appendix B Disk Profiles
The stellar disk radii that appear in Fig. 5 are generated from analysis of the distribution of stars in the simulation. Fig. 15 compares the distribution of stellar mass from the simulation with the analytic forms assumed by galform, both for disk stars (21) and for the mass of stars in the bulge, assumed to be distributed such that the projected surface density profile is given by:
Though the agreement is not comprehensive, the analytic forms are an adequate description of the global distribution of simulated stars for the majority of the systems evolution. Unsurprisingly, the simple profiles of (21) and (38) fail to describe the simulated system at early times, before an ordered galactic system has formed.
The fact that galform adheres, at all times, to the fitting forms that are really only appropriate to recent times () is clearly a valid criticism of this technique. However, as long as such forms are indeed a good description of low-redshift systems, then the only way their earlier misjudgements can significantly mislead the predictions of the model is in the early time star formation rates (Fig.7), which are based on this assumed structure.