THE AGORA HIGH-RESOLUTION GALAXY SIMULATIONS COMPARISON PROJECT

Ever since the dawn of high-performance computing, cosmological simulations have been the main theoretical tool for studying the hierarchical assembly of dark matter halos, the survival of substructure, and baryon dissipation within the cold dark matter (CDM) hierarchy, the flows of gas into and out of galaxies, and the nature of the sources responsible for the reionization, reheating, and chemical enrichment of the Universe. Purely gravitational simulations of the distribution of dark matter on large scales using different codes—e.g., Millennium I and II using Gadget (Springel et al. 2005; Boylan-Kolchin et al. 2009) and Bolshoi and BigBolshoi/MultiDark using Art (Klypin et al. 2011; Prada et al. 2012)—now produce consistent results (Springel 2012; Kuhlen et al. 2012). This is also true of collisionless "zoom-in" high-resolution simulations: Via Lactea II (using Pkdgrav-2; Diemand et al. 2008), Aquarius (using Gadget; Springel et al. 2008), and GHALO (using Pkdgrav-2; Stadel et al. 2009). To follow the formation and evolution of galaxies and clusters, however, it is necessary to model baryonic physics, dissipation, chemical enrichment, the heating and cooling of gas, the formation of stars and supermassive black holes (SMBHs), magnetic fields, non-thermal plasma processes, along with the effects of the energy outputs from these processes. A number of numerical techniques have been developed to treat gasdynamical processes in cosmological simulations, including Lagrangian smoothed particle hydrodynamics (SPH; Gingold & Monaghan 1977; Lucy 1977; Monaghan 1992) and Eulerian adaptive mesh refinement (AMR; Berger & Oliger 1984; Berger & Colella 1989). Because of the complexity of the problem, the nonlinear nature of gravitational clustering, the different assumptions made regarding the cooling and heating functions of enriched, photoionized gas, and the different implementations of crucial gas subgrid physics, it is non-trivial to validate the results of different techniques and codes even when applied to similar astrophysical problems.

The Santa Barbara Cluster Comparison project (Frenk et al. 1999) first showed the benefit of comparing hydrodynamic simulations of the same astrophysical system, a galaxy cluster, starting from the same initial conditions and using a large variety of codes. By modern standards, it was a relatively simple effort, in which the gas was assumed to be non-radiative. The spatial resolutions at the centers of simulations ranged from 5 to 400 kpc, and outputs were compared at various redshifts. The resulting simulations showed good agreement in the dark matter and gas density profiles, with a spread of about a factor of two in the predicted, resolution-dependent X-ray luminosity. Systematic differences were noted between SPH and grid-based codes including a mismatch in the central gas entropy profile, although the issue is now considered resolved (e.g., Wadsley et al. 2008; Power et al. 2013).

One of the recent comparisons of hydrodynamic galaxy simulations is the Aquila project (Scannapieco et al. 2012), in which 13 different simulations were run from the same initial conditions using various implementations of the Gadget-3 (an updated version of Gadget-2; Springel 2005) and Gasoline (Wadsley et al. 2004) SPH codes, the Arepo moving mesh code (Springel 2010), and the Ramses AMR code (Teyssier 2002). The initial conditions were those of the Aquarius halo Aq-C, with a mass at z = 0 of ∼1.6 × 10¹² M_☉, comparable to that of the Milky Way, and a relatively quiescent formation history. All groups participating in the comparison adopted their preferred implementation of radiative cooling, star formation, and feedback, and each simulation was run at two different resolutions to test convergence. None of the Aquila simulations produced a disk galaxy resembling the Milky Way. Most runs resulted in unrealistic systems with too large a stellar mass and too little cold gas, a massive bulge and a declining rotation curve. The stellar mass ranged from ∼4 × 10¹⁰ to ∼3 × 10¹¹ M_☉. Simulations with greater feedback led to lower stellar masses but usually had a hard time producing a galactic disk with a rotation velocity in agreement with the Tully–Fisher relation of late-type spirals. The star formation typically peaked at z ∼ 4 and declined thereafter, with essentially all simulations forming more than half their stars in the first ∼3 Gyr. The gaseous disk sizes were in better agreement with observations, but with too little gas mass. That the choice of numerical technique affected the Aquila results was shown, for example, by comparing Gadget-3 and Arepo runs with the similar subgrid physics implementations: the Arepo simulation produced almost twice as much stellar mass as the Gadget-3 simulation. The Aquila project shows the need to control the baryon overcooling, to prevent the early burst of star formation (e.g., Eke et al. 2000), and to promote the accretion and retention of late-accreting high-angular-momentum baryons needed to form spiral disks.

On the other hand, a new cosmological simulation of extreme dynamic range, Eris (Guedes et al. 2011), succeeded for the first time to produce a realistic massive late-type galaxy at z = 0 in which the structural properties, the mass budget in the various components (disk, bulge, halo), and the scaling relations between mass and luminosity are all consistent with a host of observational constraints. Run with the Gasoline code, Eris had 25 times better mass resolution than the typical Aquila simulation, and adopted a blast wave scheme for supernova (SN) feedback (Stinson et al. 2006) that generated galactic outflows without explicit wind particles.³⁰ Combined with a high gas density threshold for star formation, this scheme has been found to be key also in producing realistic dwarf galaxies (Governato et al. 2010). Indeed, the importance of reaching a high star formation threshold has been well studied since it was first pointed out by Kravtsov (2003). It enables energy deposition by supernovae (SNe) within small volumes and enables the development of an inhomogeneous interstellar medium (ISM) where star formation and heating by SNe occur in a clustered fashion. The resulting outflows at high redshifts reduce the baryonic content of galaxies and preferentially remove low angular momentum gas, decreasing the mass of the bulge component (Brook et al. 2011). All in all, high numerical resolution appears to be essential if simulations are to resolve the regions where stars form, and thus to succeed in producing realistic galaxies. Yet, even at Eris's resolution one barely resolves the vertical disk structure of the Milky Way since the Milky Way H i disk scale height is about 120 pc (Lockman 1984) and the H₂ scale height is about 50 pc inside the solar circle (Sanders et al. 1984; Narayan & Jog 2002).

As the above discussion should make clear, the success of cosmological galaxy formation simulations in producing realistic galaxies is a strong function of resolution (and the corresponding subgrid physics, e.g., higher star formation threshold). Numerically resolving the star-forming regions and the disk scale height is necessary because the interplay between the simulation resolution and the realism of subgrid models of star formation and feedback processes is increasingly thought of as a key to successful modeling of galaxy formation (e.g., Agertz et al. 2011; Hummels & Bryan 2012). In retrospect, this is not surprising. Stars form and deposit at least some of their feedback in the densest, coldest phase of the ISM. The characteristic property of star-forming regions is that they have extinctions high enough to block out the ultraviolet (UV) starlight that pervades most interstellar space. In the absence of UV light, the ISM undergoes a phase transition from H i to H₂ and the gas temperature drops to ∼10 K, which is likely the critical step in the onset of star formation (Krumholz et al. 2011; Glover & Clark 2012). In the Local Group, the characteristic sizes of these star-forming molecular clouds are only ∼10–100 pc. They occupy a negligible fraction of the ISM volume but contain a non-trivial fraction of its total mass: ∼30% of in the Milky Way, with lower fractions in dwarf galaxies and higher fractions in larger galaxies with denser ISMs (Blitz et al. 2007). Despite molecular clouds' high density, however, star formation within them remains surprisingly inefficient. In nearby galaxies, the observed star formation timescales in molecular gas is ∼2 Gyr (e.g., Bigiel et al. 2008; Schruba et al. 2011), and it has been known for more than 30 yr that, on average, the molecular gas forms stars at a rate of no more than ∼1% of the mass per dynamical time (Zuckerman & Evans 1974; Krumholz & Tan 2007; Krumholz et al. 2012). It is clear that any successful model for galaxy formation must include an adequate model for this critical, high-density phase, hence the need for high resolution. High resolution is also essential to avoid numerical loss of angular momentum for SPH codes that may alter the kinematics and morphology of the disk and spheroid and may lead to overly massive bulges and steep rotation curves (e.g., Kaufmann et al. 2007; Mayer et al. 2008).

The simulations in our project seek to follow the processes that regulate star formation on small scales as faithfully as possible. There are a number of feedback processes whose relative importance likely depends on the type of galaxy and the spatial scale that one is considering (see Dekel & Krumholz 2013 for a recent summary), including photoionization that heats gas to ∼10⁴ K and disperses star-forming clouds (Whitworth 1979; Matzner 2002; Dale et al. 2012), fast stellar winds that shock-heat the ISM and produce expanding bubbles (Castor et al. 1975; Weaver et al. 1977; Chevalier & Clegg 1985), the pressure of both direct starlight and dust-reprocessed radiation (Krumholz & Matzner 2009; Murray et al. 2010; Hopkins et al. 2011; Krumholz & Thompson 2012, 2013), and energy injection by Type Ia and Type II SNe. It is clear from both observations (e.g., Lopez et al. 2011) and theory (e.g., Hopkins et al. 2011; Stinson et al. 2013) that these feedback processes interact with one another in non-trivial ways. For example, ionization and momentum from radiation pressure and stellar winds act on short time and spatial scales when massive stars are formed to "clear out" the dense regions of the giant molecular cloud in which they form, and ionize and heat the surrounding neighborhood. This allows hot gas from shocked SNe ejecta—which occur several Myr later (often after the cloud is mostly destroyed)—to escape and couple to the more diffuse ISM, preventing its rapid cooling and allowing it to expand further and drive outflows.

By following these processes as directly as possible and by constraining against well-tested simulations of local systems and the ISM, one of the goals of the Assembling Galaxies Of Resolved Anatomy (AGORA) project is to lift the degeneracies between the subgrid treatments of current cosmological galaxy formation models. Indeed, the largest barrier to using today's cosmological simulations to constrain fundamental physics of cooling, shocks, turbulence, the ISM, star formation, and dark matter on sub-galactic scales is probably the unconstrained degrees of freedom in subgrid treatments of the ISM. At the same time, simulations of star formation and feedback in isolated galaxies or sub-regions of them have for the most part been run only over a narrow range in redshift, metallicity, and structural properties, without the context provided by realistic dark matter halos or baryonic inflows. The simulations to be studied in the present project provide an opportunity to explore the physics of star formation regulation in a fully cosmological context. Only by iterating between small-scale, high-resolution simulations and cosmological ones like AGORA can we hope to reach a complete theory of galaxy formation.

This motivated us to organize a new galaxy "zoom-in" simulations comparison project, with an emphasis on resolution, the physics of the ISM, feedback and galactic outflows, and initial conditions covering a range of halo masses from dwarfs to massive ellipticals. We also wanted to require more similarity in the astrophysical assumptions used in the simulations. A meeting was held at the University of California, Santa Cruz, in 2012 August to initiate this project. As of this writing, 95 individuals from 47 different institutions worldwide, using a variety of platform codes, have now agreed to participate in what has been named the AGORA project.³¹

In this paper, we describe the goals, infrastructures, techniques, and tools of the AGORA simulations comparison project. These should be of interest not only to the groups participating in AGORA but also to other groups conducting galaxy simulations, since one of the purposes of the project is to increase the level of realism in all such simulations. The numerical techniques and implementations used in this project include the SPH codes Gadget and Gasoline, and the AMR codes Art, Enzo, and Ramses (see Section 5.2 for more information). These codes will share common initial conditions (Section 2) and common astrophysics packages (Section 3), including photoionizing UV background (UVB), metal-dependent radiative cooling, metal and energy yields, stellar initial mass function (IMF), and will be systematically compared with each other and against a variety of observations using a common analysis toolkit (Section 4.2). The goals of the AGORA project are, broadly speaking, to raise the realism and predictive power of galaxy simulations and the understanding of the feedback processes that regulate galaxy "metabolism," and by doing so to solve long-standing problems in galaxy formation.

In order to achieve this goal, the AGORA project will employ simulations designed with state-of-the-art resolution. Since it is clear that the interplay between resolution and subgrid modeling of star formation and feedback is one of the key aspects of modeling galaxy formation, our choice is mandatory to make progress in this area despite the expected large computational cost. Doing this in a variety of code platforms is also essential, not only for the benefit of the groups using each code, but also to verify that the solutions are robust—i.e., that the astrophysical assumptions are responsible for any success, rather than artifacts of particular implementations. This way, the project will enable improved scientific understanding of galaxy formation, a key subject that is at last yielding to a combination of theory and computation. We plan to achieve this by sharing outputs at many redshifts, with many groups analyzing each issue using common analysis code applied to outputs from multiple codes—even timestep by timestep in some cases. Many of these intermediate timesteps along with the common initial conditions will be made available to the community.

To build the infrastructure of the AGORA project four task-oriented working groups have been established (see Table 1). These working groups ensure that the comparison of simulations is bookended by common initial conditions, common astrophysical assumptions, and common analysis tools. We also have initiated nine science-oriented working groups, each of which aims to perform original research (see Table 2; see also the AGORA project Web site for leaderships and memberships of these groups) and address basic problems in galaxy formation both theoretically and observationally. In other words, the AGORA project is not just a single set of simulations being compared, but it is a launchpad to initiate a series of science-oriented subprojects, each of which is independently designed, executed, and studied by members of the science-oriented working groups.

An example of a problem in galaxy formation that we want to address is the mechanisms that lead to galactic transformations—the processes that form galactic spheroids from disks (e.g., mergers versus disk instability), and the processes that quench star formation in galaxies (e.g., active galactic nucleus feedback versus cutoff of cold flows as halo mass increases). An important constraint that has emerged in the last few years, halo abundance matching (e.g., Behroozi et al. 2013; Moster et al. 2013), has led to a "stellar mass problem," namely, for a given halo mass, the combined mass of stellar disks and stellar bulges is too large in numerical simulations of galaxy formation relative to the expectations. Only a handful of simulations have been shown to be consistent with such constraint at z = 0 including Eris (e.g., Munshi et al. 2013), but all are seriously discrepant at higher redshift. Another difficult problem in galaxy formation is the mutual effects of baryons on dark matter, and dark matter on baryons, in accounting for the radial distribution, kinematics, and angular momentum of stars and gas in observed galaxies. It appears that compression of the central dark matter due to baryonic infall (e.g., Blumenthal et al. 1986; Gnedin et al. 2004, 2011) is an important effect in early-type galaxies, but this is apparently largely offset by other effects in late-type galaxies (e.g., Dutton et al. 2007; Trujillo-Gomez et al. 2011; Dutton et al. 2013). Thus, important astrophysical phenomena appear to cause opposite effects, which makes this a challenging problem. In order to address these processes of galaxy formation, it will be necessary to understand better the astrophysics of star formation and feedback, and the fueling and feedback from SMBHs—both of which are treated as subgrid physics in galaxy-scale simulations (e.g., Kim et al. 2011). We will tackle this challenge by carefully comparing simulations using different codes and different subgrid implementations with each other and with observations, and also by simultaneously improving the theoretical understanding of these processes, including running very high-resolution simulations on small scales.

The remainder of this paper is organized as follows. Section 2 discusses the initial conditions for the AGORA simulations, both cosmological and isolated. Two sets of cosmological initial conditions are generated using the Multi-Scale Initial Conditions (Music; Hahn & Abel 2011) code for halos with masses at z = 0 of about 10¹⁰, 10¹¹, 10¹², 10¹³ M_☉, one set with a quiescent merger history and the other set with many mergers. Additional initial conditions are generated for an isolated disk galaxy with gas fraction and structural properties characteristic of galaxies at z ∼ 1, to which the same simulation codes will be applied. Section 3 discusses the common astrophysical assumptions to be applied in all of the hydrodynamic simulations, such as the UVB, the metallicity-dependent gas cooling, and the stellar IMF and metal production assumptions. Section 4 describes strategies for running the simulations and comparing them at many redshifts, with each other and with observations. The yt analysis code (Turk et al. 2011) will be instrumental in comparing the simulation outputs, as it takes as input the outputs from all of the simulation codes being studied. The objectives of science-oriented comparison of the simulation outputs are discussed, too. Section 5 demonstrates the proof-of-concept dark-matter-only simulation for a galactic halo with a z = 0 mass of M_vir ≃ 1.7 × 10¹¹ M_☉ to field-test the pipeline of the project, including the common initial conditions and common analysis platform. The new results and ideas established in this paper are summarized in Section 6.

Common cosmological initial conditions for all simulation codes are generated using the Music code (Hahn & Abel 2011).³² Music uses a real-space convolution approach in conjunction with an adaptive multi-grid Poisson solver to generate highly accurate nested density, particle displacement, and velocity fields suitable for multi-scale "zoom-in" simulations of cosmological structure formation. For the run with two-component baryon+CDM fluid, we assume in the AGORA project that density perturbations in both fluids are equal and follow the total matter density perturbations (i.e., CDM and baryon perturbations have the same power spectrum; good approximation at late times). For particle-based codes, baryon particles are generated on a staggered lattice with respect to the CDM particles and displaced by the same displacement field as the CDM particles evaluated at the staggered positions. This strategy is to minimize two-body effects. For grid-based codes, the density perturbations are generated directly on the mesh using the local Lagrangian approximation. In both cases, the initial temperature is set to the cosmic mean baryon temperature at the starting redshift. We refer the reader to Hahn & Abel (2011) for details on the algorithms employed, but describe aspects that are particularly relevant for the AGORA project in what follows. Music generates the various fields on a range of nested levels ℓ of effective linear resolution 2^ℓ. The level covering the entire computational domain is called levelmin, ℓ_min, and the maximum level levelmax, ℓ_max.

• 1.  
  White noise generation and phase consistency. The white noise fields that source all perturbation fields are drawn reproducibly from a sequence of random number seeds {s_i}, where i ∈ [ℓ_min, ℓ_max]. The random number generator used in Music is one that comes with the GNU Scientific Library, so the white noise fields are the identical on any machine as long as they are drawn from the same seeds {s_i}. Specifying this sequence of numbers thus entirely defines the "universe" for which Music generates the initial conditions. Particular "zoom-in" regions of high resolution can be shifted or enlarged, and the resolution can be increased or decreased, without breaking consistency. This means that the Music parameter file can be distributed rather than a binary initial conditions file. Appendix A illustrates an example of such parameter files.
• 2.  
  Multi-code compatibility. Music supports initial conditions for baryons and dark matter particles for a wide range of cosmological simulation codes, many of which represented also in the AGORA project. Support for the various simulation codes and their various file formats for initial condition files is achieved through a C++ plugin mechanism. This allows adding new output formats without touching any parts of the code itself, and code specific parameters can be added transparently.
• 3.  
  Expandability. The current set-up allows expandability in various respects. (1) Support for new simulation codes can be added through plugins rather than file conversion. (2) The size and resolution of the high-resolution region can be altered consistently. (3) Future simulations focusing on larger regions or higher resolution are easily possible and will be consistent with existing simulations. (4) Cosmological models can be changed easily, and alternative perturbation transfer functions, e.g., for distinct baryon and dark matter perturbations or for warm dark matter, can be easily adopted.

Using the Music code, we have generated two sets of cosmological initial conditions for high-resolution zoom-in simulations targeted at halos with z = 0 masses of about 10¹⁰, 10¹¹, 10¹², 10¹³ M_☉, one set with a quiescent merger history (i.e., relatively few major mergers) and the other set with a violent merger history (i.e., many mergers between z = 2 and 0 for a ∼10¹⁰–10¹² M_☉ halo). Physically, these choices cover from the formation of dwarf galaxies to elliptical galaxies and to galaxy groups (see Table 3). First-order Lagrangian perturbation theory is used (default in the Music code) to initialize displacements and velocities of dark matter particles. The ΛCDM cosmological parameters we adopt are consistent with the Wilkinson Microwave Anisotropy Probe (WMAP) 7/9 results (Komatsu et al. 2011; Hinshaw et al. 2013) that includes additional cosmological data from ground-based observations of the Type Ia SNe and the baryonic acoustic oscillation: Ω_m = 0.272, Ω_Λ = 0.728, σ₈ = 0.807, n_s = 0.961, and H₀ = 70.2 km s⁻¹ Mpc⁻¹. Radiation energy density and the curvature terms are assumed to be negligible: Ω_R = Ω_k = 0. By performing a comparison study, we find that adopting the latest Planck cosmology (Planck Collaboration et al. 2013) does not noticeably change the properties of individual halos. While the present paper employs WMAP cosmology, the AGORA collaboration may later decide whether to switch to more recent cosmological parameters.

First, low-resolution dark-matter-only pathfinder simulations are performed from z = 100 to z = 0 to identify halos of appropriate merger histories. They are carried out in a (5 h⁻¹ comoving  Mpc)³ box for ∼10¹⁰ M_☉ halos, and in a (60 h⁻¹ comoving  Mpc)³ box for ∼10¹¹–10¹³ M_☉ halos. A strong isolation criterion is imposed for the quiescent set of the initial conditions to select target halos. That is, the 3R_vir radius circle of the halo being simulated must not intersect the 3R_vir radius circle of any halo with half or more of its mass at z = 0. A relaxed criterion is used for the violent set of the initial conditions: 2R_vir circle instead of 3R_vir. Then, a higher-resolution dark-matter-only simulation (e.g., particle resolution of ∼3 × 10⁵ M_☉ for ∼10¹¹–10¹³ M_☉ halos) is performed on a new initial condition re-centered on each of the target halos with nested resolution elements around it.

The highest-resolution region in this initial condition is sufficiently large to include all the structures that merge with the target galaxy or have a significant impact on its evolution. The highest-resolution region is also carefully selected so that the target halo is "contaminated" only by a minimal number of lower-resolution particles at final redshifts (z = 0 for ∼10¹⁰–10¹² M_☉ halos; z = 2 for ∼10¹³ M_☉ halos). While this region is typically a superset of a Lagrangian volume of the target halo's ∼2R_vir sphere at the final redshift, it should also be as small as possible in order to minimize the computational cost (i.e., CPU hours, memory consumption). To this end, Music supports highest-resolution particles to be placed only in a minimum bounding ellipsoid of the Lagrangian volume.³³ For particle-based codes, this determines the position of the highest-resolution region directly, while for grid-based codes, an additional refinement mask is generated that traces the high-resolution region. A high-resolution dark matter run is used to iteratively adjust the highest-resolution region by checking the contamination level inside the target halo at a final redshift. Initial conditions generated this way have been verified readable by all the participating codes in our proof-of-concept tests (see Section 5).

Stellar feedback processes are implemented quite differently in our different participating codes. Properly modeling SNe explosions or radiation from young stars remains a challenge when the target spatial resolution is 100 pc. Most, if not all, current feedback implementations are therefore highly phenomenological, and they are based on various parameters that need to be calibrated on required observational properties of simulated galaxies (e.g., star formation rate, gas and stellar fraction, H i versus stellar mass, metallicity). Moreover, most of the proposed models depend strongly on the adopted mass and spatial resolution. It is of primary importance to understand how each individual code needs to be calibrated to reproduce various observational constraints. In a comparison like the AGORA project, it is even more important to cross-calibrate stellar feedback processes of the various codes using an idealized set-up such as an isolated disk. This is precisely the goal of this second type of initial condition: we would like to model a realistic galactic disk using our various codes and their feedback recipes, varying both the feedback parameters and the mass and spatial resolutions. By doing so, subgrid star formation and feedback prescriptions in various code platforms will be tuned to provide a realistic interstellar and circumgalactic medium (CGM). We will use for that a well-defined set of observables (e.g., star formation rate, stellar and cold gas fractions, metal-enriched CGM) to identify for each code the appropriate set of feedback parameters as a function of the mass resolution. It is only after this first careful step that we will be able to move toward our final goal, namely, comparing codes in cosmological simulations.

The isolated disk galaxy initial conditions with gas fraction and structural properties characteristic of galaxies at z ∼ 1 are generated using the MakeDisk code written by Volker Springel. This code is based on solving the Jeans equations for a quasi-equilibrium multi-component halo/disk/bulge collisionless system, the particle distribution function in velocity space being assumed to be Maxwellian. Our initial conditions are quasi-equilibrium four-component systems: the dark matter halo of a circular velocity of v_c, 200 = 150 km s⁻¹ has a mass of M₂₀₀ = 1.074 × 10¹² M_☉, and follows Navarro–Frenk–White (Navarro et al. 1997) profile with corresponding concentration parameter c = 10 and spin parameter λ = 0.04. The disk follows an exponential profile (as a function of the cylindrical radius r and the vertical coordinate z) with scale length r_d = 3.432 kpc and scale height z_d = 0.1 r_d. The disk is decomposed into a stellar component of mass M_d = 4.297 × 10¹⁰ M_☉ and a gaseous component with f_gas = M_d, gas/M_d = 20%. The last ingredient is a stellar bulge that follows the Hernquist (1990) density profile with bulge-to-disk mass ratio B/D = 0.1.

We have generated three different sets of initial conditions that differ in the number of resolution elements used to describe each component (see Table 4). The low-resolution disk has 10⁵ particles for the halo, the stellar disk and the gaseous disk, and only 1.25 × 10⁴ particles for the bulge. The medium-resolution version has 10 times more particles in each component, and the high-resolution one has 100 times more elements. Note that for the gas disk, we provide particle data that can be used directly by SPH or indirectly by grid-based codes via projecting the particles into a grid. On the other hand, the default option for grid-based codes is to use the analytical density profile for the gaseous component which is just

$$\begin{equation} \rho (r,\,z)=\rho _0e^{-r/r_{\rm d}}\cdot e^{-|z|/z_{\rm d}}, \end{equation} \tag{ 1 }$$

with $\rho _0= M_{\rm d} / 4\pi r_{\rm d}^2 z_{\rm d}$. All codes will use 10⁴ K for the initial disk temperature. In grid-based codes, we also need to set the gas density, pressure, and velocity in the halo. For the halo, we recommend using zero velocity, zero metallicity, a low gas density n_H = 10⁻⁷ cm⁻³, and a high gas temperature of 10⁶ K for the halo zero velocity. This way, the total mass of gas coming from the halo is negligible compared the gas disk mass. We refer interested readers to an upcoming article in preparation for more details on the isolated disk initial conditions (e.g., the required spatial resolution for each mass resolution) and the results from the test using such initial conditions.

The rate at which diffuse gas cools radiatively determines the response of baryons to dark matter potential wells, regulates star formation, controls stellar feedback, and governs the interaction between galactic outflows and the CGM. The ejection of the nucleosynthetic products of star formation into the ISM modifies its thermal and ionization state, as radiative line transitions of carbon, oxygen, neon, and iron significantly reduces the cooling time of enriched gas in the temperature range 10–10⁷ K. The picture is further complicated by the presence of ionizing radiation, which removes electrons that would otherwise be collisionally excited and reduces the net cooling rates. Photoionization increases the relative importance of oxygen as a coolant and decreases that of carbon, helium, and especially hydrogen (Wiersma et al. 2009). Both ionizing background radiation and metal line cooling must be included for the cooling rates to be correct to within a few orders of magnitude (Tepper-García et al. 2011).

All AGORA simulations will use a standardized chemistry and cooling library, Grackle.³⁴ Grackle provides a non-equilibrium primordial chemistry network for atomic H and He (Abel et al. 1997; Anninos et al. 1997), H₂ and HD (Abel et al. 2002; Turk et al. 2009), Compton cooling off the cosmic microwave background (CMB), tabulated metal cooling and photo-heating rates calculated with the photoionization code Cloudy (Ferland et al. 2013).³⁵ GRACKLE also provides a look-up table for equilibrium cooling; depending on the problem at hand, both solvers may be used by the AGORA simulations. At each redshift, the gas is exposed to the CMB radiation and a uniform UVB and is assumed to be dust-free and optically thin. The metals are assumed to be in ionization equilibrium, so one can calculate in advance the cooling rate for a parcel of gas with a given density, temperature, and metallicity, that is photoionized by incident radiation of known spectral shape and intensity. Following Kravtsov (2003), Smith et al. (2008), Robertson & Kravtsov (2008), Wiersma et al. (2009), and Shen et al. (2010), we will use pre-computed tabulated rates from the photoionization code Cloudy at all temperatures in the range 10–10⁹ K (see Figure 1). Cloudy calculates an equilibrium solution by balancing the incident heating with the radiative cooling from a full complement of atomic and molecular transitions up to atomic number 30 (Zn). All metal cooling rates are tabulated for solar abundances as a function of total hydrogen number density, gas temperature, and redshift (as the radiation background evolves with time, see below), and are scaled linearly with metallicity. Instead of allowing Cloudy to cycle through temperatures until converging on a thermodynamic equilibrium solution, we will use the "constant temperature" command to fix the temperature externally, allowing us to utilize Cloudy's sophisticated machinery to calculate cooling rates out of thermal equilibrium. We will also deactivate the H₂ chemistry in Cloudy with the "no H₂ molecule" command since it is solved directly by the non-equilibrium chemistry solver in Grackle. Because we directly solve for the electron density and the ionization of the most abundant elements, we are able to calculate the mean molecular weight and the gas temperature.

Zoom In Zoom Out Reset image size

The default AGORA simulation will follow only the atomic gas phase, with star formation proceeding at a rate

$$\begin{equation} {d\rho _* \over dt} = {\epsilon \rho _{\rm gas} \over t_{\rm ff}} \propto \rho _{\rm gas}^{1.5} \end{equation} \tag{ 2 }$$

(i.e., locally enforcing the Schmidt law), where ρ_* and ρ_gas are the stellar and gas densities, is the star formation efficiency, and $t_{\rm ff}=\sqrt{3\pi /(32 G\rho _{\rm gas})}$ is the local free-fall time. We will also apply a density threshold below which star formation is not allowed to occur, and a non-thermal pressure floor to stabilize scales of order the smoothing length or the grid cell against gravitational collapse and avoid artificial fragmentation (Bate & Burkert 1997; Truelove et al. 1997; Robertson & Kravtsov 2008). As noted in Section 2.2, we will use non-cosmological disk simulations to tune up star formation prescription parameters for the different codes, such as the star formation density threshold, the star formation efficiency , the initial mass of star particles, and the stochasticity of star formation.

The metagalactic radiation field provides a lower limit to the intensity of the radiation to which optically thin gas may be exposed. It will be implemented in the AGORA simulations using the latest synthesis models of the evolving spectrum of the cosmic UVB by Haardt & Madau (2012). Compared to previous calculations (Haardt & Madau 1996; Faucher-Giguère et al. 2009), the new models include the following: (1) the sawtooth modulation of the background intensity from resonant line absorption in the Lyman series of cosmic hydrogen and helium, (2) the X-ray emission from the obscured and unobscured quasars that gives origin to the X-ray background, (3) an accurate treatment of the photoionization structure of absorbers that enters in the calculation of the helium continuum opacity and recombination emissivity, and (4) the UV emission from star-forming galaxies following an empirical determination of the star formation history of the Universe and detailed stellar population synthesis modeling.

The resulting UVB intensity has been shown to provide a good fit to the hydrogen-ionization rates inferred from flux decrement and proximity effect measurements, predicts that cosmological H ii (He iii) regions overlap at redshift 6.7 (2.8), and yields an optical depth to Thomson scattering that is in agreement with WMAP results. If needed, the models will be updated to include, e.g., new measurements of the mean free path of hydrogen-ionizing photons through the intergalactic medium (Rudie et al. 2013, but see also O'Meara et al. 2013).

In the AGORA simulations, each star particle represents a simple stellar population with its age, metallicity, and a Chabrier (2003) universal IMF:

$$\begin{equation} \phi (m)={dn\over dm}\propto \left\lbrace \begin{array}{@{}lc}e^{-(\log m -\log m_{\rm c})^2/2\sigma ^2}/m & (m<1\,{M_\odot })\\ m^{-2.3} & (m>1\,{M_\odot }), \end{array}\right. \end{equation} \tag{ 3 }$$

with m_c = 0.08 M_☉ and σ = 0.69. The IMF is normalized so that ∫mϕ(m)dm = 1 M_☉ between 0.1 and 100 M_☉. Star particles will inject mass and metals back into the ISM through Type II and Type Ia SNe explosions, and stellar mass loss. The time of this injection depends on stellar lifetimes in the case of Type II and on a distribution of delay times for Type Ia. The former will be determined following the Hurley et al. (2000) parameterization for stars of varying masses and metallicities.

We will follow the production of oxygen and iron; the yields of which are believed to be metallicity-independent. The masses of oxygen and iron ejected into the ISM can be converted to a total metal mass as

$$\begin{equation} M_Z = 2.09 M_{\rm O} + 1.06 M_{\rm Fe}, \end{equation} \tag{ 4 }$$

according to the solar abundances of alpha (C, N, O, Ne, Mg, Si, S) and iron (Fe, Ni) group elements of Asplund et al. (2009), as the gas cooling rate is a function of total metallicity only.

Stars with masses between 8 and 40 M_☉ explode as Type II SNe and deposit a net energy of 10⁵¹ erg into their surroundings. For the assumed IMF, the number of core collapse SNe per unit stellar mass is 0.011 M_☉⁻¹. We will use the following fitting formulae,

$$\begin{equation} M_{\rm Fe}=0.375e^{-17.94/m} \,{M_\odot } \end{equation} \tag{ 5 }$$

$$\begin{equation} M_{\rm O}=27.66e^{-51.81/m} \,{M_\odot }, \end{equation} \tag{ 6 }$$

to the total mass of ejected oxygen (including newly synthesized and initial oxygen) and iron as a function of stellar mass m (in units of M_☉) tabulated by Woosley & Heger (2007) for solar metallicity stars (see Figure 2). With the assumed IMF, the fractional masses of oxygen and iron ejected per formed stellar mass are 0.0133 and 0.0011, respectively. For comparison, the mass fractions of oxygen and iron in the Sun are 0.0057 and 0.0013 (Asplund et al. 2009). As discussed in Section 2.2, isolated disk simulations will be used by each code to tune the stellar feedback prescriptions for distributing energy and metals.

Zoom In Zoom Out Reset image size

Type Ia SNe are generally thought to be thermonuclear explosions of accreting carbon–oxygen white dwarfs in close binaries, but the nature of the mass donor star remains unknown. To determine how many SNe Ia explode at each timestep, we will adopt the most recent delay time distribution of Maoz et al. (2012). The delay times of SNe Ia are defined as the time intervals between a burst of star formation and the explosion, and follow a power-law t⁻¹ in the interval 0.1–10 Gyr. The time integrated number of SNe Ia per formed stellar mass is 0.0013 M_☉⁻¹, or about 4% of the stars formed with initial masses in the range, 3–8 M_☉, often considered for the primary stars of SN Ia progenitor systems (Maoz et al. 2012). Type Ia SNe leave no remnant, and produce

$$\begin{equation} M_{\rm Fe}=0.63\,{M_\odot }\quad M_{\rm O}=0.14\,{M_\odot } \end{equation} \tag{ 7 }$$

of iron and oxygen per event according to the carbon deflagration model W7 of Iwamoto et al. (1999).

Stars below m = 8 M_☉ return substantial fractions of their mass to the ISM as they evolve and leave behind white dwarf remnants. In all AGORA simulations, we will adopt the empirical initial-final mass relation for white dwarfs of Kalirai et al. (2008)

$$\begin{equation} w_m=(0.394+0.109m) \,{M_\odot } \end{equation} \tag{ 8 }$$

over the interval 1 M_☉ < m < 8 M_☉. We will also assume that stars with 8 M_☉ < m < m_BH = 40 M_☉ return all but a w_m = 1.4 M_☉ remnant, and stars above m_BH collapse to black holes without ejecting material into space, i.e., w_m = m. Few stars form with masses above 40 M_☉, so the impact of the latter simplifying assumption on chemical evolution is minimal. The "return fraction"—the integrated mass fraction of each generation of stars that is put back into the ISM over a Hubble time—can be written as

$$\begin{equation} R = \int _{1\,{M_\odot }}^{40\,{M_\odot }} (m-w_m)\phi (m)dm \end{equation} \tag{ 9 }$$

and is equal to 0.41 for the adopted IMF. Because the return rate is so high, stellar mass losses can prolong star formation even in systems without fresh gas inflow (e.g., Leitner & Kravtsov 2011; Voit & Donahue 2011). In practical terms, we shall implement this gas recycling mechanism by determining for each stellar particle the range of stellar masses that die during the current timestep and then calculating a returned mass fraction for this mass range using Equation (8). The metallicity of the returned gas is simply the metallicity of the star particle, i.e., we will not include metal production by intermediate mass stars.

It is worth briefly noting a few points on the common astrophysics package the AGORA project has adopted. (1) We are specifying the common astrophysics components because we want these not to be causes for inter-platform differences. It should be emphasized that we do not aim to determine "the best" models to use in galaxy simulations, nor do we attempt to undermine the freedom of choices in the numerical galaxy formation community. Any model we adopt here will be outdated soon by better theories and observations. We encourage the community to keep developing sophisticated physics and subgrid models for galaxy simulations and to investigate the problem with various methods in an independent manner. (2) Also, the AGORA common physics package is not about deciding which of the IMF, energy, or metal yields is the "correct" one. While our combination of assumptions may allow us to produce "reasonable, realistic-looking" galaxies, the sizable number of tunable and/or degenerate parameters makes it less likely to know exactly which assumptions are the "correct" ones. (3) As noted in Sections 3.2 and 3.5, we do not want to specify at this point the subgrid models for star formation and feedback which we expect to be inevitably different for each code. As Section 2.2 should make clear, we will use isolated disk initial conditions to tune the models per code, thereby constraining the inter-platform difference.

The AGORA simulations are designed and planned by the science-oriented working groups in consultation with the AGORA steering committee, while the simulations themselves will be run by experts on each participating code. The results of these runs will again be analyzed by members of several science-oriented working groups. As illustrated in Sections 1.3 and 4.3, the AGORA project is not a one-time comparison of a set of simulations, but a launchpad to initiate many subprojects, each of which is independently investigated by the science-oriented working groups. The AGORA simulations will be run and managed by the following core guidelines.

• 1.  
  Designing and running the simulations. The AGORA simulations will be designed with specific astrophysical questions in mind, so that comparing the different simulation outputs can determine whether the adopted astrophysical assumptions are responsible for any success in solving the problem in galaxy formation, rather than artifacts of particular numerical implementations. The numerical resolution is recommended to be at least as good as that of the Eris calculation in an attempt to resolve the disk scale height (Sections 1.2 and 1.3). Common initial conditions (Sections 2.1 and 2.2) and common astrophysical assumptions (Section 3) will need to be employed. Subgrid prescriptions for stellar feedback will need to be tuned to produce a realistic galaxy in an isolated set-up (Section 2.2). Resolution tests will be encouraged.
• 2.  
  Data management. Each participating code will generate large quantities of unprocessed, intermediate data, in the form of "checkpoints" describing the state of the simulation at a given time. These outputs can be used both to restart the simulation and to conduct analysis. We plan to store 200 timesteps equally spaced in expansion parameter in addition to redshift snapshots at z = 6, 3, 2, 1, 0.5, 0.2, and 0.0 at the very least. Each group is also advised to store additional outputs at slightly earlier redshift, shifted by Δz ≃ ±0.05. This extra information may be used to investigate if an inter-platform offset in time-stepping causes "timing discrepancies" for halo mergers (see Section 5.3.2 for more discussion). For many timesteps to be analyzed, central data repositories and post-processing compute time will be available at the San Diego Supercomputer Center at the University of California at San Diego, the Hyades system at the University of California at Santa Cruz, and the Data-Scope system at the John Hopkins University. Additionally, we plan to reduce the barrier to entry for the simulation data by making a subset of derived data products available through a Web interface.³⁶
• 3.  
  Public access. One of the key objectives of the AGORA project is to help interpret the massive and rapidly increasing observational data on galaxy evolution being collected with increasing angular resolution at many different wavelengths by instruments on the ground and in space. Therefore, a necessary goal of the project is providing access to both direct unprocessed data and derived data products to individuals from the broader astrophysical community. We intend to make simulation results rapidly available to the entire community, placing computational outputs on data servers in formats that enable easy comparisons with results from other simulations and with observations.

Because the simulations are being run in the same dark matter halo merger trees, it is possible to compare the cosmological evolution of target halos predicted by different codes and subgrid assumptions, sometimes halo by halo. To this end, the common analysis working group was formed (Section 1.3) to support the development of common analysis tools and to define quantitative and physically meaningful comparisons of outputs from all simulation codes. A key role will be played by the community-developed astrophysical analysis toolkit yt (Turk et al. 2011; Turk & Smith 2011; Turk 2013), which is being instrumented to natively process data from all of the simulation codes being used in the AGORA project.³⁷

• 1.  
  Science-driven analysis. yt provides a method of describing physical, rather than computational, objects inside an astrophysical simulation. For this, yt offers tools for selecting regions, applying analysis to regions, visualizing and exporting data to external analysis packages. yt allows astrophysicists to think about the physical questions, rather than the necessary computational steps to ask and answer those questions.
• 2.  
  Multi-platform analysis. In addition to the existing full support for patch-based AMR codes (Enzo), yt is starting to deliver support for analysis of octree-based AMR outputs (Art, Ramses) and particle-based outputs (Gadget, Gasoline). This way, common analysis scripts written in yt can be applied to access and investigate data from all of the simulation codes, enabling direct technology transfer between participants, ensuring reproducible scripts and results, and allowing for physically motivated questions to be asked independent of the simulation platform (see Appendix B).
• 3.  
  Open analysis. yt is freely available and open source, and it is supported by a large community of users and developers (Turk 2013), providing upstream paths for code contribution as well as detailed technical support. Any newly developed software developed in the project will be naturally available to the broader astrophysical community. Further, yt scripts and the resulting reduced data products can be shared online, enabling data analysis to be conducted by individuals regardless of their affiliation with the project.

In our proof-of-concept tests, we have demonstrated that the simulation outputs from all the participating codes of the AGORA project can be systematically analyzed and visualized in the yt platform, using unified, code-independent scripts (see Section 5). Additionally, yt can act as input for the Sunrise code (Jonsson 2006; Jonsson et al. 2010), which computes the light from simulated stellar populations (SSPs) and uses ray tracing to calculate the effects of scattering, absorption, and re-emission of this light by dust to generate mock observations and spectral energy distributions of the resulting galaxies in all wavebands.³⁸ The Sunrise outputs include realistic images of the simulated galaxies in many wavebands as FITS files, which can be compared with observed ones. It should be noted that this post-processing inevitably introduces new uncertainties including those in the SSP modeling, or in the calculation of ion number densities. Nevertheless, a detailed comparison between high-resolution simulations and the ever-increasing observational data will help constrain the numerical studies of galaxy formation.

The AGORA project will consist of a series of issue-based subprojects, each of which is studied by members of the science-oriented working groups. As shown in Table 2 of Section 1.3, these working groups intend to perform original research using multi-platform simulations and produce scientific articles for publication. Systematic and science-oriented comparisons of simulation outputs with each other and with observations are highly encouraged, not just a plain code comparison. For each science question, we will leverage the breadth of the AGORA simulations—the varied implementations of subgrid physics, hydrodynamics, and resolution—both to understand the differences between simulation outputs and to identify robust predictions. We refer the readers to the AGORA project Web site for the scientific objectives of these working groups and the project as a whole.³⁹

We design a proof-of-concept dark-matter-only simulation with a sub-L_⋆-sized galactic halo described in Section 2.1: a halo of virial mass M_vir ≃ 1.7 × 10¹¹ M_☉ at z = 0 with a quiescent merger history. For a high-resolution "zoom-in" simulation of pure dark matter, [ℓ_min, ℓ_max] = [7, 12] is selected in a (60 h⁻¹ comoving  Mpc)³ cosmological box. See Section 2.1 and Appendix A for detailed methods and parameters to generate initial conditions with the Music code. This choice of [ℓ_min, ℓ_max] corresponds to a dark matter particle resolution of 3.38 × 10⁵ M_☉ in the default highest-resolution region of 2.9 × 3.9 × 2.8  (h⁻¹ comoving  Mpc)³. In addition, using three variations of the Gadget code (Gadget-2-cfs, Gadget-3-cfs, and Gadget-3-afs; see Section 5.2.3), we have tested an initial condition in which the resolution outside the Lagrangian volume of the target halo's 2R_vir sphere is adaptively lowered (see Section 2.1 for strategies to minimize the contamination by lower-resolution particles and to set up a minimum bounding ellipsoid). The gravitational force softening length of the particle-based codes (e.g., Gadget-2/3, Gasoline, Pkdgrav-2) is set at 322 comoving pc from z = 100 to z = 9, and 322 proper pc afterward until z = 0 following Power et al. (2003). Meanwhile the finest cell size of the grid-based codes (e.g., Art-II, Enzo, Ramses) is set at 326 comoving pc, which translates into 11 additional refinement levels in a 2⁷ root grid box (Table 5). Cells of the AMR simulations are adaptively refined by factors of two in each axis on a dark matter particle overdensity of four. However, we note that the refinement algorithms used in the different AMR codes are not identical, so specifying "overdensity of four" does not fully predict the eventual refinements. Each simulation stores checkpoint outputs at multiple redshifts as described in Section 4.1 including z = 0.

We now briefly explain the participating codes in this test, focusing on the basic architecture of numerical implementations. For further information of the groups and participants using each code, we once again refer the interested readers to the AGORA project Web site. The participating codes in the future AGORA comparison studies are not limited to the ones described herein.

5.2.1. Art

5.2.2. Enzo

5.2.3. Gadget-2/3 and Their Variations

5.2.4. Gasoline and Pkdgrav-2

5.2.5. Ramses

In this section, we lay out the results of the proof-of-concept simulations. In particular, the discussion is centered on the basic analysis at the final redshift performed on the common analysis yt platform (see Appendix B for more information). More on this dark-matter-only simulation, including the detailed comparison of the halo catalogues and dark matter merger histories, will be presented in the companion paper (O. Hahn et al. 2013, in preparation; see Section 5.3.3).

5.3.1. Overall Density Structure

In Figure 3, we compile nine image panels that exhibit the results of the proof-of-concept dark-matter-only tests by nine different variations of the participating codes. Each panel displays the density-weighted projection of dark matter density in a 1 h⁻¹ Mpc box at z = 0. The overall mass distribution around the central halo shows a great similarity across all panels. The masses of the target halo are also in good agreement with one another, with $\sigma _{\rm M} / \overline{M}_{\rm vir} \sim 1.2 \%$ from the mean value, $\overline{M}_{\rm vir}$. We, however, caution that three variations of the Gadget code (Gadget-2-cfs, Gadget-3-cfs, and Gadget-3-afs; see Section 5.2.3) have employed an initial condition in which the resolution outside the Lagrangian volume of the target halo's 2R_vir sphere is adaptively lowered (see Section 5.1 for more information). At this wide field of view, large-scale tidal fields are thus inherently different depending on initial conditions and the aggressiveness of resolution choices in the lower-resolution region. Therefore, we remind the readers that particle distributions only within ∼R_vir can be most reliably compared across all nine code platforms with the best available resolution we adopted.

Zoom In Zoom Out Reset image size

For this reason, we from now on focus only on the structure in the vicinity of the central halo (R < R_vir ≃ 150 kpc). Assembled in Figure 4 are dark matter density profiles centered on the target halo of mass M_vir ≃ 1.7 × 10¹¹ M_☉ at z = 0. To make these profiles, all the dark matter particles within each radial shell are considered, including substructures and lower resolution particles, if any. As demonstrated in the bottom panel of Figure 4, all nine profiles agree very well within a fractional difference of 20% down to a radius of ∼1 kpc. The location of the maximum density is chosen to be the center of the profile; therefore, the inter-code discrepancies in the centers of profiles may explain the difference in profiles, especially within ∼1 kpc of radius. We do not find any obvious systematic difference between AMR and SPH codes, or between different gravity solvers. We again note that all the profiles in this figure are generated with a common yt script. We refer the interested readers to Appendix B to see an example script we employed for the presented analysis.

Zoom In Zoom Out Reset image size

5.3.2. Substructure Mass Distribution

Figure 5 shows the density-weighted projections of squared dark matter density of the nine different proof-of-concept runs at z = 0 in a 200 h⁻¹ kpc box. It helps to visualize where the substructures are located near the host halo within its virial radius, R_vir. Readers should note that the field of view for each panel approximately encompasses the extent of the virial radius of the target halo (R_vir ≃ 150 kpc for a M_vir ≃ 1.7 × 10¹¹ M_☉ halo). The structural differences between different code platforms in this scale are more prominent than what is observed in a wider field of view (e.g., Figure 3). The code-to-code variations at this scale could be attributed to many causes. For example, when integrated for a long time, a benignly small deviation in the density distribution at high redshift could evolve into a significant difference later and become pronounced at z = 0 especially at this highly zoomed-in scale. Because substructures on this scale are in a highly non-linear and dynamically chaotic regime, it would be unlikely to recover halo-to-halo agreements across all platforms. A relatively small timing mismatch in the numerical integration of the equations of motion could also prompt a non-negligible disparity when the runs are compared after a long integration. Indeed, the discrepancies of the effective timing of the simulations were found to be an important factor in many comparison studies, including the Santa Barbara Cluster Comparison project (Frenk et al. 1999, see also Wadsley et al. 2004 for further descriptions of the issue in the Santa Barbara comparison). These "timing discrepancy" precipitates the mismatch in the relative positions of small substructures and in the timing of substructure mergers immediately prior to z = 0 when we compare the runs.

Zoom In Zoom Out Reset image size

Another reason for code-to-code variations is the intrinsic difference in numerical methods to solve the Poisson equations for N-body dynamics. In order to quantitatively inspect such variations, substructures within 150 kpc from the center (location of the maximum density) of the target halo are identified by the Hop halo finder included in yt with an overdensity threshold δ_outer of 80 times the critical density of the Universe (Eisenstein & Hut 1998; Skory et al. 2010).⁴⁵ The resulting particle group mass functions at z = 0 are displayed in Figure 6. Only the groups containing more than 32 particles are drawn. Shown together in a dotted line is a power-law functional N(> M) = 0.01(M/M_host)⁻¹ that denotes an equal amount of mass per mass decade, to guide the reader's eye. Note that we have refrained from using the term "subhalos" to describe the particle groups identified by Hop, because the groups identified this way do not perfectly fit the typical definition of subhalos. Some of the "subhalos" within R_vir might have been linked with the host halo by the Hop algorithm.

Zoom In Zoom Out Reset image size

The close resemblances of the mass functions among the particle-based codes with tree-based gravity solvers (Gadget-2-cfs, Gadget-3, Gadget-3-cfs, Gasoline, Pkdgrav-2) and among the grid-based codes with adaptive meshes (Art-II, Enzo, Ramses) are noticeable. However, also unmistakable is the mismatch between these two breeds of codes. This phenomenon is indeed well studied and documented by many authors (e.g., O'Shea et al. 2005; Heitmann et al. 2005, 2008). They found that the AMR codes that use a multi-grid or FFT-based gravity solver achieves poorer force resolution at early times than the particle-particle-particle-mesh (P³M) or tree-PM methods in Lagrangian codes, assuming that the number of base meshes (i.e., grid cells at levelmax ℓ_max = 12 in our experiment) is roughly the number of particles, with no or little adaptive mesh at high z. Due primarily to the lack of force resolution at early redshifts, the low-mass end of the mass function tends to be suppressed for AMR codes. Consequently, it has been argued that AMR codes need more resolution in the base grid to achieve the same dark matter mass function at the low-mass end as the Lagrangian codes (e.g., O'Shea et al. 2005; Heitmann et al. 2006). Readers should note, however, the behavior of the adaptive-resolution code Gadget-3-afs, which provides results closer to the fixed-resolution codes thanks to its corrective formalism (Iannuzzi & Dolag 2011).

We emphasize that the shapes of mass functions may vary not only because of (1) the intrinsic differences in numerical techniques, but also because of (2) the inter-platform timing discrepancies discussed earlier, (3) the force and mass resolution adopted in the test, and (4) the characteristics of the halo finder. From these considerations, we argue that it would be premature, if not ill-fated, to characterize a code-to-code difference based solely on the differences in mass functions by a single halo finder at a single epoch.

5.3.3. Discussion and Future Work