THE AGORA HIGH-RESOLUTION GALAXY SIMULATIONS COMPARISON PROJECT. II. ISOLATED DISK TEST

The rate at which the gas in our galaxy radiatively cools is determined by AGORA's standard chemistry and cooling library Grackle (Bryan et al. 2014; Kim et al. 2014; Smith et al. 2016, see footnote 42). For this study, the equilibrium cooling version of Grackle is interfaced with each participating code, either via Grackle's original interface or via N. Gnedin's auxiliary API.⁵⁰ In the chosen equilibrium cooling mode, Grackle follows tabulated cooling rates pre-computed by the photoionization code Cloudy (Ferland et al. 2013).⁵¹ The pre-computed look-up table also includes metal cooling rates for solar abundances, $1\,{Z}_{\odot }$, as a function of gas number density and temperature. These metal cooling rates are then scaled linearly with metallicity which is followed in our simulations as a separate passive scalar.⁵² We also adopt metagalactic UV background radiation at z = 0 by Haardt & Madau (2012) provided by Grackle. For the difference between the chosen UV background model and previous calculations such as Haardt & Madau (1996) or Faucher-Giguère et al. (2009), we refer the readers to Section 3.3 of Kim et al. (2014).

Lastly, a non-thermal Jeans pressure floor is applied to stabilize the scales of the smoothing length (particle-based codes) or the finest cell (mesh-based codes) against unphysical collapse and to avoid artificial fragmentation due to unresolved pressure gradient (Truelove et al. 1997; Robertson & Kravtsov 2008). In practice, it is achieved by enforcing that the local Jeans length ${\lambda }_{\mathrm{Jeans}}$ be sufficiently resolved with the finest resolution elements at all times. That is,

$$\begin{eqnarray}&&{\lambda }_{\mathrm{Jeans}}={N}_{\mathrm{Jeans}}{\rm{\Delta }}x\end{eqnarray} \tag{ 2 }$$

where ${\rm{\Delta }}x=80\,\mathrm{pc}$ is the adopted spatial resolution (finest cell size or softening length; see Section 4.1) and ${N}_{\mathrm{Jeans}}=4$ is the Jeans number adopted from Truelove et al. (1997). This gives the required pressure floor value as

$$\begin{eqnarray}&&{P}_{\mathrm{Jeans}}=\displaystyle \frac{1}{\gamma \pi }{N}_{\mathrm{Jeans}}^{2}G{\rho }_{\mathrm{gas}}^{2}{\rm{\Delta }}{x}^{2}\end{eqnarray} \tag{ 3 }$$

where G is the gravitational constant, $\gamma =5/3$ is the adiabatic index, and ${\rho }_{\mathrm{gas}}$ is the gas density. Note that ${N}_{\mathrm{Jeans}}$ is not necessarily equal to the parameter controlling the pressure support in each code. For actual parameter choices for selected codes, see Appendix A. For implementations using polytropes in Art-II and Ramses, see Sections 5.2 and 5.4, respectively.

In addition to the gas physics described in the previous section, Sim-SFF incorporates subgrid models for star formation and supernova feedback. First, a parcel of gas above the threshold ${n}_{{\rm{H}},\mathrm{thres}}=10\,{\mathrm{cm}}^{-3}={\rho }_{\mathrm{gas},\mathrm{thres}}/{m}_{{\rm{H}}}$ produces stars at a rate that follows the local Schmidt law as

$$\begin{eqnarray}&&\displaystyle \frac{d{\rho }_{\star }}{{dt}}=\displaystyle \frac{{\epsilon }_{\star }{\rho }_{\mathrm{gas}}}{{t}_{\mathrm{ff}}}\end{eqnarray} \tag{ 4 }$$

where ${\rho }_{\star }$ is the stellar density, ${t}_{\mathrm{ff}}={(3\pi /(32G{\rho }_{\mathrm{gas}}))}^{1/2}$ is the local free-fall time, and ${\epsilon }_{\star }=1 \%$ is the star formation efficiency per free-fall time. We caution that the parameter ${\epsilon }_{\star }$ is not necessarily equal to the star formation efficiency parameter found in each code (e.g., for Changa and Gasoline, see Appendix B). For a new star particle to spawn, it should have at least the mass of a gas particle in the IC, ${m}_{\mathrm{gas},\mathrm{IC}}=8.593\times {10}^{4}\,{M}_{\odot }$. Note that ${n}_{{\rm{H}},\mathrm{thres}}$ adopted in this experiment is for this particular run only, and represents where the Jeans polytrope intersects with a typical $T-\rho$ equation of state in our disks (see Sections 5.2 and 5.4).⁵³

New star particles inject energy, mass, and metals back into the interstellar medium (ISM) through core-collapse (Type II) supernovae. Assuming the AGORA standard Chabrier (2003) initial mass function (IMF) and that stars with masses between 8 and $40\,\,{M}_{\odot }$ explode as Type II supernovae, one Type II supernova occurs per every $91\,{M}_{\odot }$ stellar mass formed (see Section 3.5 of Kim et al. 2014). With the AGORA recommended fitting formulae Equations (4)–(6) of Kim et al. (2014) and the assumed IMF, this single burst is found to release $2.63\,{M}_{\odot }$ of metals⁵⁴ and $14.8\,{M}_{\odot }$ of gas (including metals). Per every $91\,{M}_{\odot }$ stellar mass, these metal and mass are instantaneously deposited into its surrounding after a delay time of 5 Myr, along with a net thermal energy of 10⁵¹ erg.

We note that exact deposit schemes for energy, mass, and metals are left at each participant's discretion. We do not intend to overly specify a single common deposit scheme which will need to be inevitably different from one code to another (e.g., between mesh-based codes and particle-based codes), as we argued in Section 3.8 of Kim et al. (2014). Nevertheless, for all mesh-based codes (Art-I, Art-II, Enzo, and Ramses), the same strategy was chosen: thermal energy, mass, and metals are added to the cell where a 5 Myr old star particle sits at the time of explosion, and to this cell only. For particle-based codes (Changa, Gasoline, Gadget-3, Gear, and Gizmo), each code's deposit scheme is discussed in detail in Section 5. In future AGORA projects, we plan to calibrate different feedback schemes against observations and against one another. We refer the readers to Section 7 for more discussion on this future work.

For all codes the gas mass resolution in hydrodynamics needs to be set as close as possible to ${m}_{\mathrm{gas},\mathrm{IC}}=8.593\times {10}^{4}\,{M}_{\odot }$. Assuming that we wish to resolve a self-gravitating clump with 64 of these resolution elements, the corresponding Jeans length scale becomes

$$\begin{eqnarray}&&{\lambda }_{\mathrm{Jeans}}=2{\left(\displaystyle \frac{64{m}_{\mathrm{gas},\mathrm{IC}}}{(4\pi /3){n}_{{\rm{H}},\mathrm{thres}}}\right)}^{1/3}=348.7\,\mathrm{pc},\end{eqnarray} \tag{ 5 }$$

and therefore, from Equation (2) we choose a spatial resolution of 80 pc. This value is used as the finest cell size ${\rm{\Delta }}x$ for mesh-based codes, and as the gravitational softening length ${\epsilon }_{\mathrm{grav}}$ for particle-based codes. For all particle-based codes taking part in the present study, gravity is softened according to the cubic spline kernel (e.g., Equation (A1) of Hernquist & Katz 1989). For readers interested in the actual parameter choices, in Appendix C we examine the meanings of relevant parameters in different particle-based codes.

For particle-based codes (including Gizmo; see Section 5.9 and footnote ⁶⁷), we require that the hydrodynamical smoothing lengths for collisional particles do not drop below 20% of the gravitational softening lengths. Unlike the gravitational softening kernel, exact smoothing kernel choices differ from code to code, and are detailed for each of the particle-based codes in Section 5. We also refer the readers interested in the actual parameter choices to Appendix C again.

We recommend to mesh-based code groups that a cell be split into 8 child cells once the cell contains more mass than ${m}_{\mathrm{gas},\mathrm{IC}}=8.593\times {10}^{4}\,{M}_{\odot }$ (1 gas particle mass in the IC of particle-based codes), or 8 collisionless particles (disk/bulge star particles in the IC with ${m}_{\star ,\mathrm{IC}}=3.437\times {10}^{5}\,{M}_{\odot }$, or dark matter particles with ${m}_{\mathrm{DM}}=1.254\times {10}^{7}\,{M}_{\odot }$). This causes the grids to be refined in a fashion similar to the Lagrangian behavior of particle-based codes, and keeps the ratio of collisionless particle numbers to gas cells approximately unity on average. However, exact refinement strategies differ slightly from code to code, and are detailed for each of the mesh-based codes in Section 5.⁵⁵ We continue to refine the grids down to the resolution limit ${\rm{\Delta }}x=80$ pc (see Section 4.1) where the non-thermal pressure floor kicks in (see Section 3.1).

In Art-I, differential equations of fluid dynamics are integrated using a shock-capturing Eulerian method described in Khokhlov (1998). It uses a second-order accurate Godunov solver (Godunov 1959) that evaluates Eulerian fluxes by solving the Riemann problem at every cell interface (Colella & Glaz 1985). Left and right states of the Riemann problem are obtained by piecewise linear interpolation (van Leer 1979). In contrast to other versions of Art (see Section 5.2.1 of Kim et al. 2014), Art-I with distinctive star formation and feedback recipes (e.g., Ceverino & Klypin 2009; Ceverino et al. 2014) have been developed by A. Klypin and collaborators.

The octree-based, multi-level adaptive mesh allows users to control the grid structure at the individual cell level. For this comparison, the Art-I group uses a 128³ root grid covering a ${(1.304\mathrm{Mpc})}^{3}$ box, then achieves an ∼80 pc cell size at maximum 7 levels of refinement. The mass thresholds above which a cell is adaptively refined into an oct of 8 child cells are ${m}_{\mathrm{gas},\mathrm{IC}}=8.593\times {10}^{4}\,{M}_{\odot }$ and ${m}_{\star ,\mathrm{IC}}=3.437\times {10}^{5}\,{M}_{\odot }$ for gas and collisionless particles, respectively.⁵⁶ For the supernova feedback scheme to deposit the thermal energy adopted by all mesh-based codes, we refer the readers to Section 3.2.

Art-II solves the gravito-hydrodynamics equations using a particle-mesh + Eulerian AMR approach. Art-II features MPI parallelization for distributed memory machines, flexible time-stepping hierarchy, and a variety of unique physics modules (e.g., Gnedin & Kravtsov 2011; Agertz et al. 2013) developed by N. Gnedin, A. Kravtsov and collaborators.

For the present study, starting from a uniform 128³ root grid covering ${(1.311\mathrm{Mpc})}^{3}$, cells are refined up to 7 additional levels to reach the finest size of 80 pc. Spherical regions of 4 (6, 10) root grid cells radius around the box center are always refined to at least 3 (2, 1) additional levels relative to the root grid. The (de-)refinement procedure consists of three steps. First, cells are marked for refinement if the gas mass in the cell exceeds $0.6\,{m}_{\mathrm{gas},\mathrm{IC}}=0.6\times 8.593\times {10}^{4}\,{M}_{\odot }$, or if the cell contains two or more dark matter and/or star particles that were present in the IC. We then use a diffusion step to also mark neighboring cells for refinement and thus smooth the shape of the regions to be refined. By contrast, cells with gas masses below $0.2\,{m}_{\mathrm{gas},\mathrm{IC}}$ or without particles are marked for de-refinement provided they also satisfy a number of additional constraints. Finally, cells are refined (de-refined) by splitting them into 8 (by merging 8 child cells).

The pressure floor implemented in Art-II affects cells at the highest level of refinement by modifying the gas pressure values that enter the Riemann solver (i.e., not the actual pressure or temperature fields) with

$$\begin{eqnarray}&&{P}_{\mathrm{cell}}={\mathtt{\max }}({P}_{\mathrm{Jeans}},\,{P}_{\mathrm{gas}})\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&={\mathtt{\max }}({n}_{{\rm{H}}}{k}_{{\rm{B}}}{T}_{\mathrm{Jeans}},\,{P}_{\mathrm{gas}})\end{eqnarray} \tag{ 7 }$$

where ${P}_{\mathrm{cell}}$ is the value entering the Riemann solver, ${P}_{\mathrm{gas}}$ is the gas pressure field in the simulation, ${k}_{{\rm{B}}}$ is the Boltzmann constant, ${n}_{{\rm{H}}}={\rho }_{\mathrm{gas}}/{m}_{{\rm{H}}}$, and ${T}_{\mathrm{Jeans}}={T}_{{\rm{J}}}({n}_{{\rm{H}}}/{n}_{{\rm{H}},{\rm{J}}})$ with ${T}_{{\rm{J}}}=1800\,{\rm{K}}$ and ${n}_{{\rm{H}},{\rm{J}}}=8\,{\mathrm{cm}}^{-3}$. This polytrope choice is designed to match the common prescription Equation (3) with ${N}_{\mathrm{Jeans}}\simeq 4$. For the supernova feedback scheme to deposit the thermal energy adopted by all mesh-based codes, see Section 3.2.

Enzo is a block-structured adaptive mesh code, developed by an open-source, community-driven approach (Bryan & Norman 1997; O'Shea et al. 2004; Bryan et al. 2014).⁵⁷ Among a variety of solver choices, for this comparison the third-order accurate piecewise parabolic method is selected to reconstruct the left and right states of the Godunov problem (Colella & Woodward 1984; Bryan et al. 1995), along with a Harten–Lax–van Leer with contact (HLLC) Riemann solver (Toro et al. 1994). A maximum 30% of the required Courant–Friedrichs–Lewy (CFL) timestep is used to advance fluid elements; i.e., CFL safety factor = 0.3. In addition to solving the conservation equations for mass, momentum and energy, the equation for internal energy is also solved in parallel, and the conservative or non-conservative formulation is adaptively selected based on a local estimate of the energy truncation errors. This ensures that the gas temperature remains physical, even in highly supersonic regions.

The Enzo group uses a 64³ initial root grid covering a ${(1.311\mathrm{Mpc})}^{3}$ simulation box, then achieves 80 pc resolution with maximum 8 levels of refinement. The mass thresholds above which a cell is refined by factors of two in each axis are ${m}_{\mathrm{gas},\mathrm{IC}}\,=8.593\times {10}^{4}\,{M}_{\odot }$ and $8\,{m}_{\star ,\mathrm{IC}}=8\times 3.437\times {10}^{5}\,{M}_{\odot }$ for gas and collisionless particles, respectively (see footnote 56). The non-thermal pressure floor Equation (3) is used to modify the gas pressure inside the Riemann solver, but not to alter the actual gas energy field. For the supernova feedback scheme adopted by all mesh-based codes, we refer the readers to Section 3.2.

Ramses is an octree-based adaptive mesh code featuring an unsplit second-order accurate Monotone Upstream-centered Scheme for Conservation Laws (MUSCL) Godunov scheme for the gaseous component (Teyssier 2002).⁵⁸ For this comparison, Ramses group uses a ideal gas equation of state with $\gamma =5/3$, along with the HLLC Riemann solver (Toro et al. 1994) and the MinMod slope limiter (Roe 1986). The CFL safety factor for controlling the time step is set to 0.5. The dual energy formalism adopted in Enzo simulations (Section 5.3) is also used in Ramses runs.

For this study, starting from a uniform 128³ root grid covering ${(320\mathrm{kpc})}^{3}$, cells are refined up to 5 additional levels to achieve an ∼80 pc cell size. The refinement process works as follows. First, new refinement is triggered on a cell-by-cell basis if the baryonic mass (gas + newly formed stars) exceeds ${m}_{\mathrm{gas},\mathrm{IC}}=8.593\times {10}^{4}\,{M}_{\odot }$, or if the number of dark matter and/or star particles that are present in the IC exceeds 8.⁵⁹ We then mark additional cells by performing a mesh smoothing operation, expanding the initial area by one cell width in every direction. When new cells are created or old cells destroyed, density, momentum and internal energy are used as averaging and interpolating variables, thereby preventing a grid point with spurious temperature.

In Ramses, the gas pressure field includes the non-thermal pressure support term given by a temperature polytrope ${T}_{\mathrm{Jeans}}\ =\mu {T}_{{\rm{J}}}({n}_{{\rm{H}}}^{\prime }/{n}_{{\rm{H}},{\rm{J}}})$ with mean molecular weight μ, ${n}_{{\rm{H}}}^{\prime }={\rho }_{\mathrm{gas}}{X}_{{\rm{H}}}/{m}_{{\rm{H}}}$, ${T}_{{\rm{J}}}=1800\,{\rm{K}}$, ${n}_{{\rm{H}},{\rm{J}}}=8\,{\mathrm{cm}}^{-3}$ and ${X}_{{\rm{H}}}=0.76$.⁶⁰ As in Art-II, this polytrope approximately matches the common pressure support prescription Equation (3). Newly created star particles in Ramses have a fixed mass of ${m}_{\mathrm{gas},\mathrm{IC}}$, but they are spawned with a Poisson probability distribution whose parameters are designed to mimic the local Schmidt law, Equation (4). Lastly, for the common supernova feedback scheme adopted by all mesh-based codes, we refer the readers to Section 3.2.

Changa is a reimplimentation of Gasoline (see Section 5.6) in the Charm++ runtime system.⁶¹ Charm++ (Kale & Krishnan 1996, pp. 175–213)⁶² enables the overlap of computation and communication and provides adaptive load balancing infrastructure, allowing Changa to scale to hundreds of thousands of processor cores (Menon et al. 2015). The hydrodynamics in Changa closely follows that of Gasoline. SPH forces are calculated using the method of Ritchie & Thomas (2001), and energy is diffused using the scheme of Shen et al. (2010), both of which providing a more accurate treatment of multi-phase ISM. Timesteps are determined by the minimum of an acceleration and a CFL criterion. Furthermore, the timesteps of neighbors are kept within a factor of 2 of each other as in Saitoh & Makino (2009) in order to accurately integrate highly supersonic flows.

For this work, a kth nearest-neighbor algorithm is used to find the ${N}_{\mathrm{smooth}}=64$ nearest neighbors which are smoothed with the Wendland C4 kernel (Dehnen & Aly 2012) to determine hydrodynamic properties. Unlike conventional versions of Changa or Gasoline, the supernovae thermal energy, mass, and metals are directly distributed to the 64 neighboring gas particles.⁶³ Gas particles that are neighbors of particles that will explode as a supernova in their next timestep are put on timesteps suitable for their post-supernova thermal energy, preventing them from being on a much smaller timestep required in the CFL condition. Grackle cooling is implemented but it does not self-consistently account for the PdV work or other external sources of energy, a requirement for Changa and Gasoline's energy integration. Therefore, we split the energy integration into a half timestep of Grackle cooling, then a full timestep of PdV heating, and finally a second half timestep of cooling.

Gasoline is a massively parallel SPH code, first described in Wadsley et al. (2004), that has subsequently been updated with modern SPH features. It contains a subgrid model for turbulent mixing of metals and energy (e.g., Shen et al. 2010), a timestep limiter by Saitoh & Makino (2009) (see Section 5.5), and a geometric density estimator for SPH force expressions (see Section 2.4 of Keller et al. 2014, for a latest detailed description of the code and its performance).

For the current work, the Gasoline group uses a Wendland C4 smoothing kernel (Dehnen & Aly 2012) with ${N}_{\mathrm{smooth}}=200$ neighbors.⁶⁴ The same feedback scheme as Changa's (Section 5.5) is implemented, smoothed with the Wendland C4 kernel over 64 neighbors (not ${N}_{\mathrm{smooth}}=200$) to better match the amount of mass heated by feedback events with other particle-based codes. Gas particles that receive feedback compute their required CFL timestep at the timestep prior to receiving feedback, which helps to prevent numerical instability and overcooling. Grackle cooling is implemented by applying a half timestep of cooling, then a full timestep of external PdV heating, followed by a final half timestep of cooling, as in Changa (Section 5.5).

Gadget-3 is an updated version of Gadget-2, a cosmological tree-particle-mesh SPH code that was originally developed by V. Springel (Springel et al. 2001; Springel 2005).⁶⁵ Gadget-3 has important updates from Gadget-2, such as domain decomposition and dynamic tree reconstruction which may slightly alter the N-body dynamics. The Gadget-3 code used in this comparison is a modified version of the original Gadget-3 by K. Nagamine and his collaborators, which includes pressure-entropy formulation by Hopkins (2013), time-dependent artificial viscosity, variable smoothing lengths, among others (e.g., Choi & Nagamine 2012; Thompson et al. 2014; Aoyama et al. 2016).

For the present study, the Gadget-3 group adopts a quintic spline smoothing kernel (Morris 1996) with ${N}_{\mathrm{ngb}}=64$. The implementation of supernova feedback is based on an updated version of Todoroki (2014) that largely follows a Sedov–Taylor blast wave method outlined in Stinson et al. (2006, 2013, but not their cooling shutoff model). The exact model used in the current work is fully described in Aoyama et al. (2016) but, in brief, the implementation comprises the following steps. Every time a star particle explodes, we compute the "shock radius" based on Chevalier (1974) and McKee & Ostriker (1977), and then find the gas particles within the radius. We then inject thermal energy and metal yields into the identified gas particles within the shock radius, weighted by the SPH spline kernel. Finally, we note that the results of this version of Gadget-3 are not representative of all the Gadget-3 codes in the community, because some of the results are strongly dependent on the detailed implementations of baryonic physics, such as star formation and feedback.

Gear is a self-consistent, fully parallelized, chemo-dynamical tree SPH code (Revaz & Jablonka 2012) which is built on the publicly available Gadget-2 code (see Section 5.7; Springel 2005). The simulations reported here are run with the improvements dicussed in Revaz et al. (2016), including the pressure-entropy formulation proposed by Hopkins (2013), individual and adaptive time-stepping schemes (Durier & Dalla Vecchia 2012), artificial viscosity (Monaghan & Gingold 1983) supplemented with the Balsara switch (${f}_{{ij}}$ from Balsara 1995), and particle-based time-dependent viscosity coefficient (Rosswog et al. 2000).

For this study, the standard cubic spline smoothing kernel (Monaghan & Lattanzio 1985) in Gadget-2 is used with ${N}_{\mathrm{ngb}}=50$. The feedback energy, mass, and metal injection into the ISM is implemented following the standard SPH scheme. The implementation comprises the following steps. Every time a star particle explodes, we first find the nearest gas particles, according to the weighted number of neighbors as defined in Springel & Hernquist (2002). A desired number of neighbors ${N}_{\mathrm{ngb}}=50$ is used. Then we inject thermal energy and yields into the neighboring gas particles, weighted by the SPH spline kernel. Grackle cooling is performed after the kick step, once gas particles have eventually received supernova feedback energy and once the size of the next timestep is known. The adiabatic cooling/heating is first applied, and then the radiative one provided by Grackle.

Gizmo (Hopkins 2015) is a new mesh-free Godunov code based on discrete tracers, aimed at capturing the advantages of both Lagrangian and Eulerian techniques.⁶⁶ The numerical scheme implemented in Gizmo, initially proposed by Lanson & Vila (2008), follows the implementation of Gaburov & Nitadori (2011) and relies on the discretization of the Euler equations of hydrodynamics among a set of discrete tracers. Unlike in the moving mesh technique, where the volume is partitioned by a Voronoi tessellation, Gizmo distributes the volume fraction assigned to the tracers through a kernel function. For the current work, Gadget's standard cubic spline smoothing kernel is used with ${N}_{\mathrm{ngb}}=32$. Note, unlike SPH codes, these tracers only represent unstructured cells, sharing an "effective face" with the neighboring cells.⁶⁷ The Riemann problem is then solved across these faces using a Godunov method as in mesh-based codes, to accurately resolve shocks without artificial dissipation terms. Unlike mesh-based codes, these cells are not fixed in space and time, resulting in the scheme's Lagrangian behavior with intrinsically adaptive resolution. When the time evolution of the common face between two cells is considered, we use the second-order accurate Meshless Finite Mass method (described in Hopkins 2015). Gizmo's time-stepping scheme is fully adaptive, and closely follows Gadget-3 or Arepo (Springel 2010). It also includes a timestep limiter by Saitoh & Makino (2009) (see Section 5.5).

Gizmo's gravity solver is based on the tree algorithm inherited from Gadget-3, itself descending from Gadget-2 (see Section 5.7; Springel 2005). Gravitational softenings in Gizmo can be fixed or fully adaptive, but in the reported runs fixed softening length is used matching SPH codes. To model supernova feedback, the Gizmo simulations shown here adopts a similar feedback strategy used in Gear (Section 5.8). That is, energy, mass and metals are distributed among the neighboring gas particles/cells in a kernel-weighted fashion, but with ${N}_{\mathrm{ngb}}=32$. For star particles, timesteps are constrained to prevent supernovae from exploding in the timestep when the stars formed. Lastly, we caution that different feedback implementations using Gizmo in the literature adopt different algorithms to distribute supernova feedback energy (e.g., Hopkins et al. 2014, by the FIRE Collaboration).

We first examine the morphology of gas disks evolved in each of the codes in our experiments. In Figures 1–3, we compile nine panels that exhibit the results of the isolated disk galaxy simulations, first with radiative gas cooling but without star formation or feedback, Sim-noSF, and second with star formation and feedback, Sim-SFF, by the nine participating codes. Each panel displays the disk gas surface density in a $30\,\mathrm{kpc}$ box centered on the location of maximum gas density within 1 kpc from the center of gas mass. This centering criterion is adopted in all subsequent figures and plots. For visualizations of the particle-based codes hereafter (Figures 1–3, 14–15, 32, 34 and 35)—but not in any other analyses except these figures—yt uses an in-memory octree on which gas particles are deposited using smoothing kernels. The resolution of this octree governs the resolution of produced images. If more than eight particles are in an oct, that oct is refined into 64 child octs (i.e., yt parameters ${\mathtt{n}}\_{\mathtt{ref}}=8$ and ${\mathtt{over}}\_{\mathtt{refine}}\_{\mathtt{factor}}=2$), providing compatible or better image resolution than a typical SPH visualization. The densities are assigned to the octree in a scatter step. That is, we first calculate a particle's smoothing length, and then add the particle's density contribution to all cell centers of the octree cells that are within the particle's smoothing sphere.

We asked every code to output the state of the simulation immediately after it was initialized—so-called "0 Myr snapshot"—to allow ourselves to directly compare whether the IC generation was successful and consistent among codes. This exercise has been strenuously carried out for all the analyses items presented in Sections 6.1–6.3 and 6.7, enabling us to correct inconsistently initialized simulations early in the study. One such example, the surface density comparison of 0 Myr snapshots, is shown in Figure 1. A clear distinction in gas disk initialization between mesh-based and particle-based codes can be seen in this figure. To model the gas disk, SPH particles or Gizmo's discrete tracers are generated by drawing random numbers from the distribution function given by an analytic density profile. This by definition results in Poisson noise in the disk surface density shown here. Readers can also observe slight differences between mesh-based and particle-based codes in how the density field is represented in their calculations. By the nature of reconstructing the density from the positions of particles, the particle-based codes may smooth out the strong density contrast in the IC at the edge of the initial gas disk.

Interestingly, in Figures 2 and 3 at 500 Myr, there are other subtle differences noticeable between mesh-based and particle-based codes as well as within these sub-groups. While the peak densities and filamentary structures in the disks are very similar across all codes, it is noticeable that the mesh-based codes typically show lower densities in the inter-arm regions of the disks. The typical densities in those inter-arm regions—while not containing much mass—may differ by as much as an order of magnitude between mesh-based and particle-based codes (see also Section 6.7 and Figure 35 for a related discussion on spatial resolution). Another distinguishing aspect among the participating codes is the number of dense clumps formed. This is true in both the simulations with and without star formation and feedback (see also Section 6.4 and Figures 21 and 23 for a related comparison of newly formed stellar clumps).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 4 and 5 are the cylindrically binned gas surface density profiles for Sim-noSF and Sim-SFF, respectively. The cylindrical radius is defined as the distance from the galactic center. Raw particle fields are used for profiles of the particle-based codes, not the interpolated or smoothed fields constructed in yt. In other words, the total mass in each cylindrical annulus is divided by its area so that each gas particle contributes only to a bin in which its center falls. As noted earlier, the gas surface densities show a high degree of correspondence across all codes in both Sim-noSF and Sim-SFF. All nine profiles agree very well within a factor of a few at all radii (averaged fractional deviation for $2\lt r\lt 10$ kpc in Sim-SFF is 32.2% or 0.121 dex),⁷⁰ and can be approximated by exponential profiles at radii >1.5 kpc. Note that due to the gas consumed by star formation, the gas surface density slightly decreases from Figures 4 to 5, the latter of which could be compared with Figures 22 (newly formed stellar surface density) and 27 (star formation rate (SFR) surface density). Aside from this small decrease in density, there is little impact of supernova feedback from contrasting Sim-noSF and Sim-SFF. The inefficiency of stellar feedback is partly related to a only small increase in stellar mass in the first 500 Myr of evolution (see Sections 2 and 6.6 for more discussion).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Displayed in Figures 6 and 7 are the vertically binned gas surface density profiles for Sim-noSF and Sim-SFF, respectively. The height is defined as the absolute vertical distance from the x–y disk plane centered on the galactic center, i.e., $| {z}_{{\rm{i}}}-{z}_{\mathrm{center}}|$. Again, for particle-based codes raw particle fields are used to produce the profiles, not the interpolated or smoothed fields in yt. Note a smaller range in x-axes than in the previous figure (only one tenth of Figures 4 and 5). The surface densities begin to diverge above ∼0.6 kpc from the disk plane, but below that substantial agreement appears (averaged fractional deviation for $z\lt 0.6$ kpc in Sim-SFF is 30.4% or 0.115 dex). There is no systematic difference between mesh-based and particle-based codes, similar to what we find in radial density profiles, Figures 4 and 5.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 8 and 9 show the cylindrically binned, mass-weighted average of gas vertical heights for Sim-noSF and Sim-SFF, respectively. As defined before, the height is an absolute vertical distance from the disk plane. Thus, each line in Figures 8 and 9 represents the averaged $| {z}_{{\rm{i}}}-{z}_{\mathrm{center}}|$ as a function of cylindrical radius. Combined with Figures 6 and 7, these plots provide insight into the thicknesses of the gas disks in each of the codes. Again, no systematic difference is found between mesh-based and particle-based codes.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Here we examine the kinematics of gas disks that are evolved using each of the participating codes. First, in Figures 10 and 11 we show the gas rotation velocity curves for Sim-noSF and Sim-SFF, respectively. The curve reveals a mass-weighted average of gas rotational velocity of gas cells/particles, as a function of cylindrical radius. The cylindrical radius and rotational velocity are defined with respect to the galactic center (see Section 6.1 for our adopted centering scheme). For mesh-based codes, only the dense enough gas cells (${\rho }_{\mathrm{gas}}\gt {10}^{-25}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$) are considered in order to minimize the contribution of the gaseous halo which is present only in mesh-based codes (see Section 2 for more information about the halo gas distribution in our IC). From these two figures it is clear that there exists a very good agreement on gas kinematics in all the codes, as good as within a few percent at ∼10 kpc for both Sim-noSF and Sim-SFF (averaged fractional deviation for $2\lt r\lt 10$ kpc in Sim-SFF is 2.8% or 0.012 dex). Discrepancies in the central region (<1.5 kpc) are a result of determining the galactic center that may constantly shift its position.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figures 12 and 13 reveal the gas velocity dispersion curves for Sim-noSF and Sim-SFF, respectively. The velocity dispersion quantifies the residual velocity components of gas other than the rotational velocity found in Figures 10 and 11. In other words, each line in Figures 12 and 13 denotes a square root of a mass-weighted average of ${({v}_{{\rm{i}}}-{v}_{\mathrm{rot}}(r))}^{2}$, as a function of cylindrical radius. As in Figures 10 and 11, for mesh-based codes only dense enough cells are used to compute the dispersion. Again a good agreement is found in the velocity dispersion between all codes within a few tens of percent for both Sim-noSF and Sim-SFF (averaged fractional deviation for $2\lt r\lt 10$ kpc in Sim-SFF is 17.8% or 0.071 dex). Larger variations in the central region (<1.5 kpc) are partly due to the center determination, just like in rotation velocity curves, Figures 10 and 11. The discrepancies are also produced by vertical movement of gas in the inner disk, which is captured in mesh-based codes (Art-I, Art-II, Enzo, and Ramses) but not as well in particle-based codes (Changa, Gasoline, Gadget-3, Gear, and Gizmo) given our very particular choice of 80 pc resolution. This can be clearly seen in the bottom panels of Figures 12 and 13, in which we plot the ratio of vertical velocity dispersion (z-direction) to total velocity dispersion to illustrate the contribution by vertical movement of gas. From these figures, we find that a significant portion of gas velocity dispersion in the inner disk measured in mesh-based codes are driven by vertical movement, but not in particle-based codes (see Section 6.7 and Figure 35 for a related discussion on spatial resolution).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figures 14 and 15, we compile 9 panels of the density-square-weighted gas temperature projections for Sim-noSF and Sim-SFF, respectively.⁷¹ Each panel is for the central $30\,\mathrm{kpc}$ box (see Section 6.1 for our adopted centering scheme). Readers should note that, by design, only mesh-based codes initially include a gaseous halo with low density and high temperature (${T}_{\mathrm{gas}}={10}^{6}\,{\rm{K}};$ colored dark red), but not the particle-based codes. All of these codes show similar features without star formation and feedback in Figure 14. We continue a related discussion on Sim-noSF using Figures 16 and 18 later in this section.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

But slight differences in supernova feedback implementation—even within a common guideline—may give rise to different features in the temperature map; Figure 15. For example, small hot bubbles are visible in Changa and Gear, but not in Gasoline or Gadget-3 in which the common supernova feedback schemes are implemented slightly differently.⁷² We also note the chimneys of hot gas departing from the disk in Enzo and Ramses, clearly visible in the edge-on maps of Figure 15. This gas ejecta from the disk as a result of supernova feedback is not as evident in Art-I, Art-II, or particle-based codes, although some hot bubbles are seen in Changa and Gear. We caution that the spatial resolution employed in this paper is only 80 pc. It is not as high as the resolutions in some of the modern zoom-in cosmological simulations, and may not be enough for particle-based codes to resolve the chimney-like structure above and below the disk (see Section 6.7 and Figure 35 for a related discussion on spatial resolution). We refer the readers to a continued discussion on Sim-SFF using Figures 17 and 19 later in this section.⁷³

Zoom In Zoom Out Reset image size

To better understand what we find in Figures 14 and 15, we show the two-dimensional probability distribution functions (PDFs) of gas density and temperature in Figures 16 and 17 for Sim-noSF and Sim-SFF, respectively. We consider gas within 15 kpc from the galactic center. Colors represent the total gas mass in each two-dimensional bin. As explained in Section 6.1, raw particle fields are used for the PDFs of particle-based codes, not the smoothed fields constructed by yt. Note that a gaseous halo, represented by low-density, high-temperature gas in the upper left corner of each panel, is designed to exist only in mesh-based codes, but is absent in particle-based codes (SPH codes and Gizmo; see Section 2). It is therefore not the intended scope of this paper to compare this hot halo or circumgalactic medium between codes. To guide the eye, we plot the mean temperature in each density bin from Changa's Sim-noSF run with a thick dashed line in each panel (in both Figures 16 and 17; in the range of $[{10}^{-26},{10}^{-21}]\,{\rm{g}}\,{\mathrm{cm}}^{-3}$). Changa's mean profile is close to the "mean of means" of these nine codes, thus it helps to compare the PDFs' relative positions between codes. The thin dotted diagonal lines denote the slope of constant pressure process, and the thin dotted–dashed diagonal lines that of constant entropy process.

Overall, all codes exhibit similar behaviors in Figure 16 when without star formation and feedback, just like the broad similarity observed in Figure 14. A clear branch of gas is visible in all codes extending toward higher-density, lower-temperature, owing to the common treatment of cooling by the Grackle library (see also Figure 18, and a related discussion later in this section). But, as noted in Section 3.1, due to varying Grackle versions participating codes are interfaced with (Grackle v2.1 in Changa, Gasoline, Gadget-3 and Gizmo, versus Grackle v2.0 or below in Art-I, Art-II, Enzo, Ramses and Gear), the cooling rates differ slightly from code to code even at the same density, temperature, and metal fraction. See Section 2 (footnote ⁴⁸) or Section 3.1 (footnote ⁵²) for more information. Any remaining discrepancy is attributable to the difference in how each code sub-cycles its cooling module in the hydrodynamics calculation.

Zoom In Zoom Out Reset image size

Now we compare the density–temperature PDFs when star formation and supernova feedback are included in Figure 17. All codes successfully lower the fraction of low-temperature, high-density gas by forming stars and then injecting their feedback energy (see also Figures 18–20, and a related qualitative discussion later in this section). However, notable differences exist between codes as to how gas reacts to the supernova feedback. For example, in Changa and Gear some gas is leaving the aforementioned high-density branch toward higher temperature (up to 10⁶ K) due to supernova feedback, but not in Gasoline, Gadget-3 and Gizmo. The high-temperature, high-density gas seen in Changa and Gear is associated with the small hot bubbles discussed in Figure 15. As explained earlier, these discrepancies in particle-based codes are attributed to different numerical implementations of the common feedback physics. Also noticeable is the hot gas being ejected from the disk as a result of feedback particularly in Enzo and Ramses, seen as a broader distribution of hot gas in Figure 17 compared to Figure 16.

For a more qualitative comparison of the ISM thermal structure, shown in Figures 18 and 19 are the gas mass distributions along the density axis (density PDF) for Sim-noSF and Sim-SFF, respectively, simply derived from Figures 16 and 17 above. Note again that a gaseous halo, represented by low-density tails toward the left side of these plots, is by design included only in the mesh-based codes, but not in the particle-based codes. In Figure 18, when without star formation and feedback all codes show similar distributions within a factor of a few difference in the density range $[{10}^{-25},{10}^{-22}]\,{\rm{g}}\,{\mathrm{cm}}^{-3}$. A notable deviation is that three particle-based codes, Gadget-3, Gear and Gizmo, hold more mass at density above ${10}^{-22}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ than the rest of the codes do. However, in Figure 19, now in a more realistic setup including star formation and feedback, no clear systematic difference exists between mesh-based and particle-based codes. While the agreement in the range $[{10}^{-25},{10}^{-22}]\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ is again within a factor of a few (averaged fractional deviation in ${10}^{-25}\lt \rho \lt {10}^{-22}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ is 28.6% or 0.109 dex), the codes diverge from one another up to more than an order of magnitude at density above ${10}^{-22}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$, translating into differences in clumping properties (Figure 23) and SFRs (Figure 26).

Zoom In Zoom Out Reset image size

Then, Figure 20 plots the code-by-code mass change in each density bin from Figures 18 to 19. This plot aims to measure the effect of star formation and supernova feedback by subtracting the density probability distribution of Sim-noSF from that of Sim-SFF. As noted in our discussion of Figure 17, Sim-SFF lowers the fraction of high-density gas by forming stars at above ${n}_{{\rm{H}},\mathrm{thres}}=10\,{\mathrm{cm}}^{-3}={\rho }_{\mathrm{gas},\mathrm{thres}}/{m}_{{\rm{H}}}$ (see Section 3.2) and then injecting their feedback energy. This impact is more evident in the bottom panel of Figure 20, where we show the sign-preserving logarithm, symlog(), of the mass change, making smaller changes more discernible.⁷⁴ Without exception, gas masses at density above ${\rho }_{\mathrm{gas},\mathrm{thres}}=1.67\times {10}^{-23}\,{\rm{g}}\,{\mathrm{cm}}^{-3}$ (thick dashed line) are reduced by star formation and feedback, even if inefficiently. The gas is either consumed by star formation or redistributed by thermal supernova feedback to a less dense region below ${\rho }_{\mathrm{gas},\mathrm{thres}}$ away from star-forming sites.

Zoom In Zoom Out Reset image size

In this section, we study the morphology of stellar disks formed in Sim-SFF. Figure 21 shows the distributions of newly formed star particles in Sim-SFF for each of the 9 codes. Only the newly formed star particles are drawn, not the disk or bulge stars that were present in the IC. Star particles present in the IC are excluded from all "stellar" particle analyses hereafter. Each frame is centered on the galactic center, defined in Section 6.1 as the location of peak gas density within 1 kpc from the center of gas mass. This center almost always coincides approximately with the center of the most massive stellar clump. Colors represent the total newly formed stellar masses in each two-dimensional bin. The bottom rows also highlights clumps of newly formed stars (more discussion on stellar clumps later in this section).

Zoom In Zoom Out Reset image size

Figure 22 depicts surface densities of newly formed star particles—excluding star particles present in the ICs—for Sim-SFF, calculated in cylindrically symmetric radial bins. While there are differences up to a factor of a few among some codes (e.g., between Gadget-3 and Gasoline, due to their different rates of galaxy-wide star formation shown in Figure 26; averaged fractional deviation for $2\lt r\lt 10$ kpc is 53.9% or 0.187 dex), all the lines can be well fit by an exponential disk profile at radii >1.5 kpc. Occasional fluctuations visible in the profiles (e.g., at ∼6.5 kpc in Gizmo) are due to dense stellar clumps located at these particular radii.

Zoom In Zoom Out Reset image size

In order to compare the distribution of newly formed star particles and the level of disk fragmentation between different codes, we identify clumps in the distribution of newly formed star particles using the friends-of-friends (FOF) algorithm (Efstathiou et al. 1985).⁷⁵ We only consider clumps with newly formed stellar masses above $2.6\times {10}^{6}\,{M}_{\odot }$ (equivalent to 30 times the mass of star particles present in the IC, $30\,{m}_{\star ,\mathrm{IC}}$). The most massive clump found by FOF is excluded since it is always associated with the stellar bulge at the galactic center, but we do not explicitly remove gravitationally unbound clumps (which may have overestimated the number of clumps). Identified clumps are marked with circles in the bottom row of Figure 21, where the radius of each circle indicates the clump's virial radius. We also show the cumulative mass functions of newly formed stellar clumps in Figure 23. It is worth noting that very similar clump distributions and mass functions are discovered when we identify clumps using the HOP algorithm (Eisenstein & Hut 1998) in lieu of FOF. The general trends for clumps discussed below are largely independent of the clump finder.

Zoom In Zoom Out Reset image size

In all codes, the majority of newly formed stellar clumps have masses below ${10}^{7.5}\,{M}_{\odot }$, and there is a relatively sharp decline in the number of clumps toward higher masses.⁷⁶ From these figures, it is also clear that nearly formed stellar clumps are the most prevalent in the Gizmo run, and less so in the Ramses run than in other codes. The relatively large number of stellar clumps in Gizmo is not a transient feature, but consistently observed across snapshots until 1 Gyr. This is related to the fact that Gizmo produces the most clumpy gas disk among the codes even with star formation and feedback in Sim-SFF (see Figure 3). While preserving all the common elements in comparison (such as pressure floor or resolution; described in Sections 3 and 4), the Gizmo group has carried out extensive tests with other simulation parameters to check what most dictates the level of fragmentation (e.g., smoothing kernel, ${N}_{\mathrm{ngb}}$, Grackle version, slope limiter, dual energy formalism; but always within conventional norms), finding that the different parameter choices do not qualitatively alter the mass function. However, we note that increasing the Jeans pressure floor by as little as 100%—well within the uncertainty in the geometry prefactor of Jeans pressure support equation, Equation (3)—and reverting to Grackle v2.0—where metal cooling rates are slightly lower compared to v2.1 because of the different solar metallicity definition—entirely removes the discrepancy (see the black dashed line in Figure 23 labeled Gizmo-ps2). A possible explanation for the stronger fragmentation in Gizmo is that the "effective" gravitational resolution in Gizmo is slightly higher compared to other codes (as observed by Few et al. 2016), resulting from a combination of different choices in their implementation (e.g., slope limiter, gradient estimator, density estimator).⁷⁷

Following the analysis shown in Section 6.2 for the gas disk kinematics, here we study the kinematics of stellar disks formed by each code in Sim-SFF. In Figure 24 we show the rotation velocity curves for newly formed star particles in Sim-SFF. As in gas rotation velocity curves of Figures 11, each line represents a mass-weighted average of stellar rotational velocities as a function of cylindrical radius. As in Section 6.4, only the newly formed star particles are considered for these profiles, not the disk or bulge stars in the IC. Just as in the gas rotation curves, the stellar rotation velocities show a high degree of similarity among the different codes, as good as within a few percent at certain radii (averaged fractional deviation for $2\lt r\lt 10$ kpc is 2.5% or 0.011 dex). Disagreements seen at radius >12 kpc for some codes such as Art-II and Gasoline are attributed to a small number statistics near the edge of a stellar disk.

Zoom In Zoom Out Reset image size

In Figure 25 we analyze the velocity dispersion curves for newly formed star particles in Sim-SFF. As in gas velocity dispersion curves of Figures 13, this plot quantifies the residual velocity components other than the rotational velocity computed in Figure 24. Once again, a good agreement is found that all codes lie within a few tens of percent from one another at radii <10 kpc (averaged fractional deviation for $2\lt r\lt 10$ kpc is 11.2% or 0.046 dex). When compared with Figures 12 and 13, the agreement is particularly better in the central region (<1.5 kpc). Note that when we plot the ratio of vertical velocity dispersion to total velocity dispersion in the bottom panel of Figure 25, the systematic discrepancy found in Figure 13 between mesh-based and particle-based codes at radii <1.5 kpc no longer exists. This confirms our assertion in Section 6.2 that the said discrepancy in Figure 13 is due to vertical gas movement captured only in mesh-based codes.

Zoom In Zoom Out Reset image size

In this section, we compare the SFRs of different codes in Sim-SFF and check whether we reproduce the observed relation between gas surface density and SFR surface density. In the following discussion, SFRs are time-averaged over the past 20 Myr and derived based on the ages of newly formed star particles in the 500 Myr snapshots.

Figure 26 displays the evolution of the galaxy-wide SFR of each run by 500 Myr. For most codes, SFRs increase at early times, reach a maximum at 200–300 Myr after the start of the simulation, and then plateaus at later times. The SFRs evolve smoothly without evidence for strong bursts. The SFRs are within a factor of ∼3 from one another at all times (averaged fractional deviation for $50\lt t\lt 500$ Myr is 32.8% or 0.123 dex), and for most codes the values settle at 2–4 $\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ over most of the simulated time. Some differences are noteworthy. For example, Gasoline and Enzo predict somewhat lower SFRs, especially at intermediate times, while the SFR for Gizmo never plateaus or begins to decline, but reaches a maximum of $\sim 6\,\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ at 500 Myr. Gadget-3 produces the most stellar mass in this time period, but its SFR does not further grow after ∼300 Myr. The total stellar mass formed in 500 Myr ranges from $0.8\times {10}^{9}\,{M}_{\odot }$ in Gasoline to $2.4\times {10}^{9}\,{M}_{\odot }$ in Gadget-3. Gizmo's efficient and clumpy star formation is discussed in detail in Section 6.4 and Figure 23. We note again that increasing the Jeans pressure floor in Gizmo by a factor of 2—well within the uncertainty in the geometry prefactor of Jeans pressure support equation, Equation (3)—and reverting to Grackle v2.0 entirely removes Gizmo's discrepancy (see the black dashed line in Figure 26 labeled Gizmo-ps2; compare with Figure 23). In addition to the systematic tests the Gizmo group performed, the Gadget-3 group has also tried a run with a slight variation in the treatment of supernova feedback that would increase the number of gas particles inside the "shock radius" (Section 5.7). We find that this slight variation indeed produces significantly lower SFR for Gadget-3, closer to Enzo's value (results not shown here). Our tests strongly suggest that the SFR evolution is highly sensitive to the details of the numerical implementation of the common subgrid physics, including pressure floor and feedback prescriptions.

Zoom In Zoom Out Reset image size

Now, to better understand the differences in SFR, we plot in Figure 27 the SFR surface densities as functions of cylindrical distance from the galactic center. Again, SFRs are estimated using the newly formed star particles that are younger than 20 Myr old. The agreement between the codes is generally encouraging, especially outside the central 0.5 kpc. We note, however, that the SFR within 0.5 kpc constitutes a large fraction of the total galactic SFR. Enzo and Gasoline have significantly lower SFRs in the central region of the galaxy, thus explaining the difference seen in Figure 26. In contrast, much of the excess star formation in Gadget-3 and Gizmo takes place at large radii and is likely related to the formation of larger numbers of stellar clumps, as seen in Figures 21 and 23. For example, a ∼6.5 kpc peak in Gizmo's SFR profile in Figure 27 coincides with a peak at the same radius in its stellar surface density profile in Figure 22.

Zoom In Zoom Out Reset image size

When measured on galactic scales (∼kpc), the gas surface density and SFR surface density are tightly correlated (e.g., Schmidt 1959; Kennicutt 1989, 1998; Kennicutt et al. 2007; Bigiel et al. 2008). Consequently, the relation between these two observables, the Kennicutt–Schmidt (KS) relation, is frequently used to calibrate the modeling of star formation in galaxy-scale simulations. In Figure 28 we show the KS relation for Sim-SFF, where gas and SFR surface densities are measured within cylindrical annuli, as computed in Figures 5 and 27, respectively. Only annuli with nonzero averaged SFR surface densities are considered. The thick black dashed line denotes a best observational fit by Equation (8) in Kennicutt et al. (2007), $\mathrm{log}\,{{\rm{\Sigma }}}_{\mathrm{SFR}}=1.37\,\mathrm{log}\,{{\rm{\Sigma }}}_{\mathrm{gas}}-3.78$, for a spatially resolved patches in M51a. The blue hatched contour marks the observed sub-kpc patches in nearby galaxies taken from Figure 8 of Bigiel et al. (2008), where their hydrogen surface density is multiplied by 1.36 to account for helium to match the total gas density in our simulations (see their Section 2.3.1).⁷⁸

Zoom In Zoom Out Reset image size

As Figure 28 reveals, all participating codes predict a KS relation that agrees well with one another within a factor of a few, and with observed nearby disk galaxies in Bigiel et al. (2008). In particular, by combining star formation and thermal supernova feedback, most codes match both the normalization of the observed relation and the characteristic "threshold" value of the gas surface density ($\sim 10\,\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$) below which star formation becomes less efficient. However, it should be noted that our simulations do not include multi-phase gas physics that explicitly models the transition between atomic and molecular hydrogen at $\sim 10\,\,{M}_{\odot }\,{\mathrm{pc}}^{-2}$ (e.g., Krumholz et al. 2009). More investigation may be needed to check how the apparent change in slope seen here is affected by our choice of subgrid physics (e.g., star formation efficiency ${\epsilon }_{\star }$, or thermal feedback energy budget; see Section 3.2). Figure 28 also highlights that there are some differences between the codes. SFR surface densities of Gasoline and Enzo lie slightly below the other codes at a given gas surface density, while Gadget-3 and Gizmo show higher SFRs than the rest of the codes. These differences are generally in line with what was observed in global SFR of Figure 26.

In order to better match the observational technique by Bigiel et al. (2008), one may consider using sub-kpc patches to generate the KS relation rather than cylindrical annuli. In Figures 29 and 30 we present mock observations of gas surface densities and SFR surface densities for Sim-SFF at 500 Myr. For mesh-based codes in Figure 29, their panels in Figure 3 are degraded to 750 pc resolution. For gas particles for particle-based codes in Figure 29 and young star particles of age less than 20 Myr in Figure 30, we use the cloud-in-cell (CIC) scheme to deposit particle masses on to a uniform two-dimensional grid with 750 pc resolution.⁷⁹ This resolution matches Bigiel et al. (2008)'s reported working resolution. The dynamic range of the color axis are set to 3 orders of magnitude in both Figures 29 and 30, in order to help readers to see if the gas depletion times are similar among pixels.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Then we identify all 750 × 750 pc patches with nonzero SFR surface density, and plot them in Figure 31 on the same KS plane as Figure 28. As an example, all such patches found in the Changa run are shown as gray triangles. For all other codes, only 80% percentile contours are drawn. Again, all participating codes reproduce the slope and normalization of the observed KS relation well. But slight differences in SFR surface densities exist. Contours for Gasoline and Enzo lie below the pack, while Gadget-3 and Gizmo's contours sit above all other codes at all gas surface densities. These findings are consistent with what we see in Figures 26 and 28.

Zoom In Zoom Out Reset image size

In Figure 32 we present the projections of density-square-weighted gas metal fraction for Sim-SFF with star formation and feedback. The metal fraction we show here is simply the ratio of metal density to total gas density, and it is not in units of ${Z}_{\odot }$, in order to minimize any confusion caused by Grackle 2.0 versus 2.1 implementations (see footnotes ⁴⁸ and ⁵²). The color axis spans from 0.01 to 0.04 in a linear scale. The edge-on views of mesh-based codes, particularly Enzo and Ramses, show high metallicity filaments flowing out of the disk, carrying metals into the embedding halo (see also Figure 15). However, as noted in Section 2, a gaseous halo exists only in mesh-based codes, but not in particle-based codes (SPH codes and Gizmo). Therefore, it is not the intended scope of this paper to compare the metal content of the halo which, by design, is captured only in mesh-based codes. In the face-on views within the disk we find qualitatively similar results across codes, with denser gas (corresponding to star-forming regions) tending to be more metal-enriched. Significant differences in the morphology of metal distribution exist, however. Differences in numerical implementations of the stellar feedback model are responsible for such discrepancies (see Section 6.3 for more discussion).⁷³

Zoom In Zoom Out Reset image size

Figure 33 shows the mass-weighted averages of gas metal fraction in Sim-SFF on a two-dimensional density–temperature plane for gas within 15 kpc of the galactic center. Raw particle fields are used for particle-based codes, not the interpolated or smoothed fields constructed in yt. Note again that a gaseous halo represented by low-density, high-temperature gas in the upper left corner of each panel exists only in mesh-based codes. For particle-based codes, high-density, low-temperature gas has higher metallicity because of its correspondingly higher SFR and thus metal enrichment. For mesh-based codes, on the other hand, high metal fraction is found in low-density, high-temperature gas as well, which is contaminated by hot metal-enriched materials dispersed by supernova feedback.⁸⁰ These two observations are related to how well metals get mixed in each code. Readers may have noticed the difference between mesh-based and particle-based codes already in Figure 32 by focusing on mixing of the metals. With neither the halo gas nor a sophisticated metal mixing scheme in place, metal-enriched gas in particle-based codes tends to stay near the dense star-forming sites that provided the metals.

Zoom In Zoom Out Reset image size

In Figure 34 is the density-weighted average of gas elevations from the x–y disk plane for Sim-SFF with star formation and feedback. That is, averages of $({z}_{{\rm{i}}}-{z}_{\mathrm{center}})$ such that a positive (negative) value indicates the gas along that line of sight is located above (below) the disk plane on average. For particle-based codes, we use the reconstructed density field from yt's in-memory octree on which gas particles are deposited using smoothing kernels (see Section 6.1). This figure helps to visualize and estimate the warping of the gas disk (e.g., Levine et al. 2006). All participating codes produce largely flat disks, with vertical offsets less than ±1 kpc. Yet, it is also true that all codes show some levels of coherent warping or flaring along the disk plane. This strongly suggests that all these codes are able to resolve vertical instabilities.

Zoom In Zoom Out Reset image size

Finally, Figure 35 compares the sizes of spatial resolution elements each code imposes on the galactic disk of Sim-SFF at 500 Myr. For mesh-based codes, this is a projection of gas cell sizes along different lines of sight, weighted by (cell volume)⁻² so that the maximal resolution element along that line of sight could stand out. For particle-based codes, this is a projection of gas particle sizes, defined as ${({m}_{\mathrm{gas}}/{\rho }_{\mathrm{gas}})}^{1/3}$, smoothed on to yt's octree (the same octree yt used in other edge-on/face-on visualizations; see Section 6.1 for more information on the octree), weighted by (particle volume)⁻², defined as ${({m}_{\mathrm{gas}}/{\rho }_{\mathrm{gas}})}^{-2}$. The color axis spans from 10 to 10³ pc in a logarithmic scale, with highest resolution shown in dark blue. The green color in mesh-based codes marks the finest mesh size permitted, 80 kpc, while the blue color in the spirals and clumps of particle-based codes can be associated with the minimum smoothing length permitted, $0.2\times 80$ pc. Because of its Lagrangian nature, particle-based codes best demonstrate their strengths in dense clumps and spirals. Meanwhile, with its flexible refinement strategy, mesh-based codes may best utilize their strengths in the contact regions with high density contrast, such as above and below the disk. For example, in the edge-on views of Figure 35, the green-colored high resolution region covering the disk is thicker in mesh-based codes than in particle-based codes.⁸¹

Zoom In Zoom Out Reset image size