Gas Cooling in Hydrodynamic Simulations with An Exact Time Integration Scheme
Abstract
We implement and test the exact time integration method proposed by Townsend (2009) for gas cooling in cosmological hydrodynamic simulations. The errors using this time integrator for the internal energy are limited by the resolution of the cooling tables and are insensitive to the size of the timestep, improving accuracy relative to explicit or implicit schemes when the cooling time is short. We compare results with different time integrators for gas cooling in cosmological hydrodynamic simulations. We find that the temperature of the gas in filaments before accreting into dark matter halos to form stars, obtained with the exact cooling integration, lies close to the equilibrium where radiative cooling balances heating from the UV background. For comparison, the gas temperature without the exact integrator shows substantial deviations from the equilibrium relation. Galaxy stellar masses with the exact cooling technique agree reasonably well, but are systematically lower than the results obtained by the other integration schemes, reducing the need for feedback to suppress star formation. Our implementation of the exact cooling technique is provided and can be easily incorporated into any hydrodynamic code.
keywords:
galaxies: formation – evolution – methods: numerical1 Introduction
Energy loss via radiative cooling in hot temperature plasmas plays a crucial role in the formation and evolution of galaxies (Blumenthal et al., 1984). Hydrodynamic simulations are routinely used to study galaxy formation, including radiative cooling, star formation, and other relevant processes. Radiative cooling itself together with gravity in cosmology represents a complex problem that cannot be addressed analytically (e.g. Birnboim & Dekel, 2003; Kereš et al., 2005; Nelson et al., 2016).
Recently, many studies have focused on the accuracies of hydrodynamic solvers (Springel, 2010; Bauer & Springel, 2012; Hopkins, 2015; Zhu et al., 2015) and the differences between different codes in galaxy formation studies (e.g., Sijacki et al. 2012; Kereš et al. 2012; Vogelsberger et al. 2012; Torrey et al. 2012; Nelson et al. 2013; Zhu & Li 2016). While stateoftheart codes such as Arepo and Gizmo are secondorder accurate in time (and space), the cooling step is often only firstorder. It is thus also timely to review the performance of the cooling method now being used in current hydrodynamic simulations.
Typically, radiative cooling in cosmological hydrodynamic codes is either handled with an explicit scheme
(1) 
or an implicit one
(2) 
depending on whether the starting time or the ending time , with the timestep, is used to evaluate the cooling rates. Equations (1, 2) operate on the thermal energy per unit mass for gas element , a common thermodynamic state used in numerical simulations. Note that in these two Equations and the rest of the paper, we have normalized the cooling/heating rate with a factor to be consistent with the units in the cooling table.
The explicit scheme is straightforward to implement. However, stability requires that the timestep be controlled carefully. The integration of the specific internal energy in Grackle (Smith et al., 2017) is close to a firstorder explicit method. The implicit method for radiative cooling is common among SPH codes (Hernquist & Katz, 1989), in particular Gadget (Springel et al., 2001), as well as in AMR codes such as RAMSES (Teyssier, 2002). While this approach is stable, it entails inaccuracies due to the finite timestep and errors arising from multiple zero points for certain rooting finding methods (Townsend, 2009).
To control the accuracy or to maintain the stability of explicit techniques, the internal energy update in each subcycle is subject to empirical constraints such as the “10 percent rule” (Teyssier, 2015; Smith et al., 2017). For explicit schemes, subcycling is mandatory due to its unstable behavior with respect to the timestep. However, substantial errors result once the timestep is comparable to the cooling time even when subcycling is applied. This can be seen in Figure 1, which shows the evolution of the temperature in the constant density cooling test of Smith et al. (2017) using a density of . We use the cooling routines available in Grackle and Gadget, respectively. We vary the timestep to see how the temperature of the gas evolves from as a function of time. The cooling rate is calculated for metalfree gas exposed to the UV background of Haardt & Madau (2012) at redshift . For a timestep comparable to , errors appear in the predicted temperature using the Gadget cooling method. The accuracy is greatly improved with Grackle due to subcycling. Nevertheless, fluctuations in the gas temperature around are present, which is still noticeable for ^{1}^{1}1 In the case of net heating, we find that the impact of timestep size on the evolution of the gas temperature is also present for both implicit and explicit schemes..
We also repeat a test in Townsend (2009) (their Figure 1 and 2), which is shown in the right panel with the predicted gas temperature as a function of a single timestep. This test represents the situation where the timestep size can span a large range in hydrodynamics simulations. The implicit method cools the gas too quickly once the timestep size is larger than 8.5 Myr. While the gas temperature is overall stable, it is slightly above the result shown in the left panel. Thanks to subcycling, gas temperature with Grackle agrees well with the result with very fine timesteps. However, the oscillations in gas temperature appear more irregular than that in the middle panel.
For the explicit and implicit methods, the evolution of internal energy (gas temperature) is treated as an ODE system. A new integration scheme proposed by Townsend (2009) is based on the following observation. For gas cooling, unlike a usual ODE system, we already have the full knowledge of all the future evolution, which is simply expressed by the cooling curves. Thus, one can avoid the barrier presented by the stiff source terms and analytically or numerically integrate the change of internal energy along the cooling curves instead. Using this approach, under the assumption that the cooling curves depend on the temperature data points (power law or linear), one can simply and uniquely map an array of internal energies to a time array.
In this paper, we implement and extend the method of Townsend (2009) to gas cooling in cosmological hydrodynamic simulations, including cooling functions for gas exposed to a UV background. We describe a set of simple algebraic equations to be solved in Section 2. In Section 3, we test the validity of the algorithm and confirm that the accuracy of the exact time integration does not rely on the timestep size. We then apply this cooling method to a cosmological hydrodynamic simulation with gas cooling and star formation, and compare the results using the Gadget cooling method and Grackle. We discuss our findings in Section 4 and provide conclusions in Section 5. Our implementation of the exact cooling technique is publicly available.
2 Methods
We discuss an implementation of the exact integration scheme with the cooling/heating curves computed with the photoionization code CLOUDY (Ferland et al., 2013). A redshiftdependent UV ionizing background is supplied to cloudy to model the heating rate due to young stars and AGNs. The essential part of “exactintegration” is to integrate the possibly stiff cooling/heating source term numerically by assuming that these rates vary linearly^{2}^{2}2For the assumption of powerlaw dependence, we refer the reader to Townsend (2009). between the temperature grid points. Thus, numerical errors associated with this integration scheme entirely arise from the above assumption, so that they are of the same order as the linear interpolation errors. Such an approach, compared with the popular implicit time integration method, shows greatly improved accuracy. In light of the recent developments of hydrodynamic solvers in computational galaxy formation and evolution studies, we here consider whether or not errors due to radiative cooling are important remaining sources of error that should be addressed when designing progressively more sophisticated models of star formation and feedback.
For each particle or cell , the approach first locates its position in redshiftdensity grids and then bilinearly interpolates the cooling / heating rates along all the temperature data points ranging from . Note that contributions for both primordial gas composition (H+He) and metalenriched gas can be calculated. In our implementation, we have combined the mean molecular weight and gas temperature into the specific internal energy through
(3) 
such that the interpolations are carried out between the specific internal energy and
(4) 
a dimensionless temporal evolution function introduced by Townsend (2009) in their Equation (24). This function measures the total time it takes to cool the gas normalized by a reference cooling time.
We denote the upper internal energy in the cooling table as and the cooling rate at that temperature as . We use at for the adopted cooling table, but it can be chosen at any temperature higher than the current state. The computation of can be performed in a piecewise linear fashion on all the grid data points. In each segment , we assume the cooling rate varies linearly from to as
(5) 
The change in in such a segment can be analytically calculated as
(6) 
or
(7) 
In Equation (6) or Equation (7), we have defined and as
(8) 
A monatomic curve of can be now summed up cumulatively using
(9) 
which starts from a reference internal energy to current internal energy .
To compute the internal energy for gas element at the next step , we find the expected change in as:
(10) 
Following Equation (10), we identify the segment in that contains . The final step is to analytically invert Equation (6) or Equation (7) to find the predicted internal energy at the next time step. The above approach can be trivially extended to a net heating regime. One just starts from the lowest temperature as in the cooling table to compute and other quantities.
3 Results
The above equations are straightforward to implement and we provide an example C program^{3}^{3}3https://goo.gl/UqKv68, which can be incorporated into codes such as Gadget or Gizmo.
3.1 Validation of the exact integration scheme for gas cooling
Before we use this method in actual simulations, we first apply it to the constant density cooling test similar to Smith et al. (2017), which examines the accuracy of the solver over a period of time. A gas cell is at a number density of and an initial temperature. The gas cooling rate is calculated for gas exposed to the UV background of Haardt & Madau (2012) at redshift .
In the top panel of Figure 2, we show the gas temperature and the mean molecular weight as a function of a single timestep for a starting temperature of . The gas temperature is expected to rise up to due to the heating effect of the UV background, which completes at 1 Myr. The predicted temperature of the metalfree gas case is very similar to gas.
The middle panel shows the gas temperature and for a starting temperature of K, which is the same as that in Smith et al. (2017). The gas temperature and for gas agree perfectly with Figure 3 in Smith et al. (2017). In the case of net cooling, metalfree gas cools more slowly than , which is also subject to metal cooling. Note that cooling from to K is a very abrupt phase, much shorter than what it takes from to K due to heating.
The above tests show that the exact integration algorithm and our implementation are both correct in the regimes of net heating and net cooling. Moreover, the temperature is insensitive to the timestep size, which contrasts with the behavior of the explicit and implicit integrations. The predicted gas temperature shows a sharp drop for , which is very difficult to capture accurately with either explicit or implicit methods if the timestep is comparable to the cooling time.
The lower panel of Figure 2 shows a zoomedin view of the gas temperature and in the cooling test in the middle panel around the time of the sharp temperature drop. The predicted values of temperature and mean molecular weight accurately reflect the features in the cooling curves between and K (a local minimum between two peaks) within a short duration of 0.2 Myr.
3.2 Cosmological hydrodynamic simulations with different time integration schemes for gas cooling
scheme name  Gadget  Grackle  exact integration 

equation  Eq. (2)  Eq. (1)  Eq. (10) 
order of accuracy  first order  first order  exact 
subcycling  no  yes  no 
Having established that we have a correct implementation of the exact integration, we now use cosmological hydrodynamic simulations to evaluate the differences between the exact integration of gas cooling with an implicit method as in the Gadget code and the method in Grackle. We generate initial conditions in CDM using the music code (Hahn & Abel, 2011). A uniform box of comoving Mpc on each side contains a total number of gas particles and dark matter particles. Hydrodynamics is evolved using the movingfinitemass (MFM) method in the Gizmo code. Physical processes include radiative cooling and star formation. In order to isolate the effect of radiative cooling due to time integration, we do not include any “kinetic feedback”. Only the thermal pressure in the ISM (Springel & Hernquist, 2003) is applied.
The first simulation uses the exactintegration from the previous section with the cooling curves from “CloudyData_UVB=FG2011.h5” in Grackle. The cooling curves are based on CLOUDY calculations with the UV background of FaucherGiguère et al. (2010a). Our second simulation adopts a first order backward Euler method implemented in Gadget. To solve Eq. (2), Gadget uses a bisection root finding method. The last simulation is carried out with Grackle in its tabulated cooling mode also with identical cooling curves as in the other two simulations. Similar to Springel et al. (2001), we apply a restriction such that a particle is only allowed to lose at most 50% of its internal energy in any timestep for all three integration schemes. The numerical parameters controlling the gravity calculation, hydrodynamics, and star formation are identical in the three simulations so that any differences are caused by the time integration schemes. The differences between the methods are listed in Table 1.
In the top panels of Figure 3, we show the gas density distribution in the simulation box at redshift simulated with the cooling method in Gadget, Grackle, and exact integration scheme in the entire simulation box. Gas temperature is colorcoded such that cold gas appears in blue while hot gas is in yellow. All three simulations produce almost identical structures in the distribution of gas density and temperature. The bottom panels further show the zoomedin view of the gas distribution in the most massive dark matter halo. Again, differences in the distribution of gas density and temperature in the three simulations are not apparent.
3.2.1 Gas density–temperature phase diagram
It is not surprising to find that the impact of time integration for gas cooling appears to be slight due to the fact that gas cooling times in lowdensity regions and in shocked heated regions (galaxy groups) are both very long. It is imperative to consider a close inspection of the gas temperature in the regime where cooling is efficient.
To do this, it is a common practice in the field to examine the gas densitytemperature diagram. To calculate gas temperature from the specific internal energy, we use Eq. (3) with the mean molecular weight in the Grackle cooling table as a function of density, the specific internal energy, and redshift. This step ensures that the conversion from gas internal energy to temperature is accurate^{4}^{4}4Since the cooling term is treated with an operator splitting approach, the update of the mean molecular weight occurs between the two kick operations. There is a slight mismatch between the gas internal energy and the mean molecular weight in the snapshots. while the differences in the gas density–temperature phase diagrams are solely determined by the time integration.
In Figure 4, we show the gas density–temperature phase diagrams. In these plots, we also emphasize the region where heating due to the UV background dominates in red and the region of net cooling is in blue. The boundary between net cooling and net heating is shown in a solid gray line, which is an equilibrium where the heating rate equals the cooling rate. Since the intensity and the spectrum of the UV background evolve as a function of redshift, the equilibrium line also varies accordingly. Between and , the equilibrium between heating and cooling drops from for gas at density to for .
Overall, the gas phase diagrams obtained with the three integration schemes agree well in most of the parameter space. However, the temperature for gas in the density range is slightly higher in the exact time integration than the implicit integration and Grackle, which contains many gas particles with a temperature below K. Physically, this part of the phase diagram corresponds to gas falling into dark matter haloes along large scale filaments.
The gas temperature obtained with the cooling method in Gadget shows clear departures from equilibrium at and in . In particular, the distribution of the gas temperature is skewed below K for gas at densities around . Interestingly, the gas temperature follows the equilibrium relation very well at . For the simulation using Grackle cooling, the gas temperature is consistently lower than the equilibrium line in the density range .
For comparison, the gas temperature obtained with the exact integration closely follows the equilibrium relation at all three redshifts, as shown in Figure 4. Moreover, the distribution of gas particles with temperature above and below the equilibrium line is symmetrical.
The coupling between the gas dynamics and gas cooling in the hydrodynamic simulations adds more complexity in the interpretation of the departure from the equilibrium temperature with Gadget and Grackle in the phase diagram compared to the constant density cooling test in Figure 1. As the gas in the cosmic web enters into the gravitational potential of galaxies, it cools towards the equilibrium temperature but gas density continuously increases which moves the equilibrium temperature lower. The increasingly higher density could introduce a bias which amplifies the numerical error towards a lower gas temperature as the gas density can even reach a higher value.
In the case of Grackle, another bias towards a lower temperature is also present. The cooling rate and the heating rate are not symmetric around the equilibrium temperature. This bias is present in Figure 2 as it takes much shorter time to cool the gas from to K than heat it from to K. Once the temperature oscillates around the equilibrium, one would expect to find more gas particles in net heating than in net cooling.
3.2.2 – relation
While gas cooling and star formation proceed almost the same in the three simulations, some quantitative differences are to be expected in the stellar mass using different time integration schemes. Overall, we find there are actually minor differences in the three star formation histories which are reflected in the total stellar mass produced in the simulation box. The stellar mass obtained with Grackle cooling is slightly higher than that obtained with the cooling method in Gadget. The exact integration scheme produces the least stellar mass among the three. Overall, the differences in total stellar mass are small, at the level of tens of percents. Nevertheless, this is an encouraging and welcoming feature of the exact integration scheme as it suffers least from the overcooling problem.
We look further at the stellar mass in individual dark matter haloes. In Figure 5 we compare the stellar mass–dark matter mass relations at with different time integration schemes. In this plot, the scattered data points show the dark matter halo mass and the stellar mass for each FOF halo and the solid lines are locally weighted regression results.
We use the stellar mass from the exact integration as a baseline for comparison. The differences in stellar mass is shown in the lower part of Figure 5 in terms of . The stellar mass for a given dark matter halo mass in the three simulations agrees to within dex. Nevertheless, the differences in stellar mass between the exact approach and the other integrators is systematic for dark matter haloes more massive than , with the exact technique yielding the lowest masses. Since the peak of star formation efficiencies occurs at halos, the exact integration scheme could at least alleviate the quantitive need for supernovae and/or AGN feedback required in the current hydrodynamic simulations (Vogelsberger et al., 2013; Pillepich et al., 2017; Weinberger et al., 2017).
3.2.3 Is the cooling time well resolved in our simulations?
Inevitably, our hydrodynamic simulations are constrained by both mass and force resolutions. Another resolution effect, not discussed often in the literature, is worth some discussion here. Even with a set of welldefined cooling curves and accurate hydrodynamic solvers available today, the large range of cooling times in a typical cosmological simulation itself presents a significant numerical challenge. Short gas cooling times on the order of yr in Wiersma et al. (2009) indicates that there could be certain gas phases not wellresolved in current cosmological simulations.
In numerical simulations, many codes do not apply any restrictions on timesteps due to radiative cooling. We use the system timestep from the simulations, which is mostly determined by the CFL condition (Saitoh & Makino, 2010). In our tests, we have applied a restriction that each particle could at most lose 50% its thermal energy in a single step. We would expect this constraint to become fully redundant once the shortest cooling time is fully resolved.
Figure 6 shows the distribution of the ratio between the cooling time and the timestep for each particle at with the exact time integration. We use a simple mapping from the ratios to colors of gas particles as shown in the colorbar. A very long cooling time with respect to is shown in green while very short cooling times are shown in blue. The distribution of the ratio between the two time timescales spans almost ten orders of magnitude.
For the shockheated region, the low density gas with temperature above K, the cooling time is much longer than the timestep. Gas particles near the equilibrium where the cooling rate equals the heating rate also show long cooling times, which is expected since the net cooling rate approaches zero. This is why the phasespace diagrams in the three simulations appear similar to each other.
The equilibrium where the cooling and heating rates are equal shows up in Figure 6 as a prominent green line at and . For gas particles in the filaments ( and ), the gas temperature lies around the equilibrium in the exact integration scheme. On the other hand, in this phase does not exceed as the green line terminates around . This is caused by more dynamical states in the filaments where it is difficult for gas particles to lie perfectly on the equilibrium as they fall into the gravitational potential and interact with ambient gas. Above and below the equilibrium line, there is a population of gas particles with . This short timescale for gas cooling translates into a very large single timestep, which is challenging for both Gadget and Grackle cooling methods to maintain their accuracy. This is also why there are differences in the gas density–temperature phase diagrams in the filaments. Note that temperature around the equilibrium line with the exact cooling scheme could still contain substantial numerical errors due to the fact that none of the three integration schemes resolve the cooling time well. Ultrahigh resolution simulations (e.g., Nelson et al., 2016) are essentially needed to better understand the gas thermal and dynamical states in this phase.
4 Discussion
The exact time integration scheme we have studied in this work can be easily extended to more sophisticated applications to include selfshielding (Rahmati et al., 2013), or a local radiation field (Vogelsberger et al., 2013; Gnedin & Hollon, 2012) or molecular hydrogen () cooling, which requires one to generate cooling tables using photoionization codes such as CLOUDY. On the other hand, this method is difficult to apply to nonequilibrium cooling.
scheme name  Gadget  Grackle  exact integration 
average walltime ^{5}^{5}5The average time is taken directly from the measurement in raw walltime statistics. Due to the differences between the actual nodes carrying out the simulations, a more accurate comparison would be the fraction of each procedure within a single timestep.  178 s  75 s  87 s 
gravity  20.8 s  18.9 s  20.1 s 
(tree+PM)  12%  25%  23% 
34.2 s  29.7 s  33.9 s  
hydrodynamics  19%  40%  39% 
domain  8.1 s  9.3 s  9.6 s 
decomposition  5%  12%  11% 
cooling  107.6 s  7.1 s  9.4 s 
and star formation  60%  9%  11% 
The exact time integration is efficient when compared to the implicit integration which invokes multiple calls to interpolate gas cooling rates. In Table 2, we compare the walltime statistics of the three integration schemes used in our simulations. The walltime spent by gravity, hydrodynamics, and domain decomposition in the three simulations are almost the same. The major difference is in gas cooling: the exact integration consumes one tenth as much as the implicit scheme while the latter is significantly slowed by multiple table lookups and interpolations, which we have not optimized, during its root finding step. In theory, the cost of the exact integration scales as , the product of the number of gas particles and the number of temperature grid points .
This time integration scheme might also help to resolve gas temperature for studies of cooling emission (Haiman et al., 2000; Fardal et al., 2001; FaucherGiguère et al., 2010b). Because gas at temperatures around K in the current simulations is not wellresolved due to the short cooling time and large intrinsic variations in the cooling timescale, as we showed in the previous section, large uncertainties in the cooling emission could be present as the Lyman emissivity changes exponentially around K (FaucherGiguère et al., 2010a).
While the differences in the stellar mass produced by the three time integration schemes for gas cooling are relatively small ( dex), the galaxy formation model we used in our simulation does not resolve the ISM dynamics once the gas is above the threshold density (just slightly above ) for star formation. As we see from the cooling timescale distribution in Figure 6, gas cooling for most of the particles in our simulations is wellresolved except for gas more dense than . The limited range in gas densities in our simulations where there are large errors in gas cooling could be the reason why only relatively small differences in the stellar mass are found. One might expect that larger differences due to different time integration schemes could be present, as shown in the right panel of Figure 1, for simulations aiming at explicitly resolving the ISM physics.
5 Conclusions
In this paper, we have extended the method of exact time integration by Townsend (2009) to use a redshiftdependent cooling table with a UV background. We have applied this method to a cosmological hydrodynamic simulation and compared its performance with other time integration schemes. The impact of numerical errors with the current integration in gas cooling on galaxy formation appears slight but is present nevertheless. Based on our findings, we argue that the exact time integration for gas cooling is attractive because

It is able to map the specific internal energy to the timestep while it is also insensitive to the timestep. Errors are only from interpolation of the cooling table.

It is easy to implement and efficient to execute.

It is also very flexible in that it can be adapted to different cooling tables which incorporate different UV backgrounds, metaldependent cooling, and other physical processes such as including local radiation fields.

It could in principle moderate the quantitive need for supernovae and/or AGN feedback in the current models.
Acknowledgements
We thank the referee for his/her thorough review with very helpful comments on earlier drafts of the manuscript. This work is supported by NSF grants AST0965694, AST1009867, and AST1412719. We acknowledge the Institute For CyberScience at The Pennsylvania State University for providing computational resources and services that have contributed to the research results reported in this paper. The Institute for Gravitation and the Cosmos is supported by the Eberly College of Science and the Office of the Senior Vice President for Research at the Pennsylvania State University. Some of the computations in this paper were run on the Odyssey cluster supported by the FAS Division of Science, Research Computing Group at Harvard University.
References
 Bauer & Springel (2012) Bauer A., Springel V., 2012, MNRAS, 423, 2558
 Birnboim & Dekel (2003) Birnboim Y., Dekel A., 2003, MNRAS, 345, 349
 Blumenthal et al. (1984) Blumenthal G. R., Faber S. M., Primack J. R., Rees M. J., 1984, Nature, 311, 517
 Fardal et al. (2001) Fardal M. A., Katz N., Gardner J. P., Hernquist L., Weinberg D. H., Davé R., 2001, ApJ, 562, 605
 FaucherGiguère et al. (2010a) FaucherGiguère C.A., Kereš D., Dijkstra M., Hernquist L., Zaldarriaga M., 2010a, ApJ, 725, 633
 FaucherGiguère et al. (2010b) FaucherGiguère C.A., Kereš D., Dijkstra M., Hernquist L., Zaldarriaga M., 2010b, ApJ, 725, 633
 Ferland et al. (2013) Ferland G. J., et al., 2013, Rev. Mex. Astron. Astrofis., 49, 137
 Gnedin & Hollon (2012) Gnedin N. Y., Hollon N., 2012, ApJS, 202, 13
 Haardt & Madau (2012) Haardt F., Madau P., 2012, ApJ, 746, 125
 Hahn & Abel (2011) Hahn O., Abel T., 2011, MNRAS, 415, 2101
 Haiman et al. (2000) Haiman Z., Spaans M., Quataert E., 2000, ApJ, 537, L5
 Hernquist & Katz (1989) Hernquist L., Katz N., 1989, ApJS, 70, 419
 Hopkins (2015) Hopkins P. F., 2015, MNRAS, 450, 53
 Kereš et al. (2005) Kereš D., Katz N., Weinberg D. H., Davé R., 2005, MNRAS, 363, 2
 Kereš et al. (2012) Kereš D., Vogelsberger M., Sijacki D., Springel V., Hernquist L., 2012, MNRAS, 425, 2027
 Nelson et al. (2013) Nelson D., Vogelsberger M., Genel S., Sijacki D., Kereš D., Springel V., Hernquist L., 2013, MNRAS, 429, 3353
 Nelson et al. (2016) Nelson D., Genel S., Pillepich A., Vogelsberger M., Springel V., Hernquist L., 2016, MNRAS, 460, 2881
 Pillepich et al. (2017) Pillepich A., et al., 2017, preprint, (arXiv:1703.02970)
 Rahmati et al. (2013) Rahmati A., Pawlik A. H., Raičevic̀ M., Schaye J., 2013, MNRAS, 430, 2427
 Saitoh & Makino (2010) Saitoh T. R., Makino J., 2010, PASJ, 62, 301
 Sijacki et al. (2012) Sijacki D., Vogelsberger M., Kereš D., Springel V., Hernquist L., 2012, MNRAS, 424, 2999
 Smith et al. (2017) Smith B. D., et al., 2017, MNRAS, 466, 2217
 Springel (2010) Springel V., 2010, MNRAS, 401, 791
 Springel & Hernquist (2003) Springel V., Hernquist L., 2003, MNRAS, 339, 289
 Springel et al. (2001) Springel V., Yoshida N., White S. D. M., 2001, New Astron., 6, 79
 Teyssier (2002) Teyssier R., 2002, A&A, 385, 337
 Teyssier (2015) Teyssier R., 2015, ARA&A, 53, 325
 Torrey et al. (2012) Torrey P., Vogelsberger M., Sijacki D., Springel V., Hernquist L., 2012, MNRAS, 427, 2224
 Townsend (2009) Townsend R. H. D., 2009, ApJS, 181, 391
 Vogelsberger et al. (2012) Vogelsberger M., Sijacki D., Kereš D., Springel V., Hernquist L., 2012, MNRAS, 425, 3024
 Vogelsberger et al. (2013) Vogelsberger M., Genel S., Sijacki D., Torrey P., Springel V., Hernquist L., 2013, MNRAS, 436, 3031
 Weinberger et al. (2017) Weinberger R., et al., 2017, MNRAS, 465, 3291
 Wiersma et al. (2009) Wiersma R. P. C., Schaye J., Smith B. D., 2009, MNRAS, 393, 99
 Zhu & Li (2016) Zhu Q., Li Y., 2016, ApJ, 831, 52
 Zhu et al. (2015) Zhu Q., Hernquist L., Li Y., 2015, ApJ, 800, 6