Gas Cooling in Hydrodynamic Simulations with An Exact Time Integration Scheme

TL;DR

This paper implements and tests an exact time integration method for gas cooling in cosmological simulations, demonstrating improved accuracy and more realistic gas temperature and star formation results compared to traditional schemes.

Contribution

The authors adapt Townsend's exact cooling integration method for cosmological hydrodynamics, showing it enhances accuracy and reduces the need for feedback in galaxy formation models.

Findings

Gas temperature near equilibrium with radiative cooling and UV heating is accurately captured.

The exact integrator produces systematically lower stellar masses, reducing feedback dependence.

Errors are limited by cooling table resolution and are insensitive to timestep size.

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.

Figures

Tables

Table 1. Table 1: A comparison of the three time integration schemes for gas used in our cosmological hydrodynamic simulations.

scheme name	Gadget	Grackle	exact integration
equation	Eq. (2)	Eq. (1)	Eq. (10)
order of accuracy	first order	first order	exact
sub-cycling	no	yes	no

Table 2. Table 2: A comparison of wall-time spent by three time integration schemes.

scheme name	Gadget	Grackle	exact integration
average wall-time 555The average time is taken directly from the measurement in raw wall-time 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%

Equations

$${u}_i(t+\Delta t)=u_i(t)+\Delta t\cdot \Lambda (u_i(t))$$

$${u}_i(t+\Delta t)=u_i(t)+\Delta t\cdot \Lambda (u_i(t+\Delta t)),$$

$$u=\frac{kT}{(\gamma -1)\mu m_{\mathrm{H}}},$$

$$Y(u)=\frac{{\Lambda }_{\mathrm{ref}}}{u_{\mathrm{ref}}}{\int }_u^{u_{\mathrm{ref}}}\frac{\mathrm{d}u}{\Lambda (u)},$$

$$\Lambda (u)=\frac{{\Lambda }_2-{\Lambda }_1}{u_2-u_1}(u-u_1)+{\Lambda }_1,\text{if}{\Lambda }_2\ne {\Lambda }_1.$$

$$\delta Y=\frac{{\Lambda }_{\mathrm{ref}}}{u_{\mathrm{ref}}}A\log (\frac{B+u_2}{B+u_1}),\text{if}{\Lambda }_2\ne {\Lambda }_1,$$

$$\delta Y=\frac{{\Lambda }_{\mathrm{ref}}}{u_{\mathrm{ref}}}\frac{u_2-u_1}{{\Lambda }_1},\text{if}{\Lambda }_2={\Lambda }_1.$$

$$\begin{gathered} A\equiv \frac{u_2-u_1}{{\Lambda }_2-{\Lambda }_1}, \\ B\equiv \frac{u_2{\Lambda }_1-u_1{\Lambda }_2}{{\Lambda }_2-{\Lambda }_1}. \end{gathered}$$

$$Y(u)=\sum _{u_{\mathrm{ref}}}^u\delta Y,$$

$$\Delta Y=\frac{u_i(t)}{{\Lambda }_{u_i}(t)}\frac{{\Lambda }_{u_{\mathrm{ref}}}}{u_{\mathrm{ref}}}\frac{\Delta t}{t_{\mathrm{cool}}}$$

Full text

††pagerange: Gas Cooling in Hydrodynamic Simulations with An Exact Time Integration Scheme–References††pubyear: 2016

1 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 state-of-the-art codes such as Arepo and Gizmo are second-order accurate in time (and space), the cooling step is often only first-order. 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

[TABLE]

or an implicit one

[TABLE]

depending on whether the starting time $t$ or the ending time $t+\Delta t$, with $\Delta t$ the timestep, is used to evaluate the cooling rates. Equations (1, 2) operate on the thermal energy per unit mass $u_{i}$ for gas element $i$, 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 $\Lambda$ with a factor ${n_{\rm{H},i}^{2}/\rho_{i}}$ 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 first-order 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 sub-cycle is subject to empirical constraints such as the “10 percent rule” (Teyssier, 2015; Smith et al., 2017). For explicit schemes, sub-cycling 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 $t_{\rm cool}$ even when sub-cycling 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 $0.1\,\rm{cm}^{-3}$. 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 $10^{6}\,\rm{K}$ as a function of time. The cooling rate is calculated for metal-free gas exposed to the UV background of Haardt & Madau (2012) at redshift $z=0$. For a timestep comparable to $t_{\rm cool}$, errors appear in the predicted temperature using the Gadget cooling method. The accuracy is greatly improved with Grackle due to sub-cycling. Nevertheless, fluctuations in the gas temperature around $10^{4}\,\rm{K}$ are present, which is still noticeable for $\Delta t=0.01\,t_{\rm cool}$111 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 sub-cycling, 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 redshift-dependent UV ionizing background is supplied to cloudy to model the heating rate due to young stars and AGNs. The essential part of “exact-integration” is to integrate the possibly stiff cooling/heating source term numerically by assuming that these rates vary linearly222For the assumption of power-law 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 $i$, the approach first locates its position in redshift-density grids $[z,\log(n_{\rm H})]$ and then bi-linearly interpolates the cooling / heating rates along all the temperature data points ranging from $[10,10^{9}]\,\rm{K}$. Note that contributions for both primordial gas composition (H+He) and metal-enriched gas can be calculated. In our implementation, we have combined the mean molecular weight $\mu$ and gas temperature $T$ into the specific internal energy $u$ through

[TABLE]

such that the interpolations are carried out between the specific internal energy $u$ and

[TABLE]

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 $u_{\rm ref}$ and the cooling rate at that temperature as $\Lambda_{\rm ref}$. We use $u_{\rm ref}$ at $10^{9}\,\rm{K}$ for the adopted cooling table, but it can be chosen at any temperature higher than the current state. The computation of $Y(u)$ can be performed in a piece-wise linear fashion on all the grid data points. In each segment $[u_{1},u_{2}]$, we assume the cooling rate varies linearly from $\Lambda_{1}$ to $\Lambda_{2}$ as

[TABLE]

The change in $\delta Y$ in such a segment can be analytically calculated as

[TABLE]

or

[TABLE]

In Equation (6) or Equation (7), we have defined $A$ and $B$ as

[TABLE]

A monatomic curve of $Y(u)$ can be now summed up cumulatively using

[TABLE]

which starts from a reference internal energy $u_{\rm ref}$ to current internal energy $u$.

To compute the internal energy for gas element $i$ at the next step $t+\Delta t$ , we find the expected change in $Y(u_{i}(t))$ as:

[TABLE]

Following Equation (10), we identify the segment in $Y$ that contains $Y(u_{i})+\Delta Y$. 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 $u_{\rm ref}$ in the cooling table to compute $Y(u)$ and other quantities.

3 Results

The above equations are straightforward to implement and we provide an example C program333https://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 $0.1\,\rm{cm}^{-3}$ and an initial temperature. The gas cooling rate is calculated for $0.1Z_{\odot}$ gas exposed to the UV background of Haardt & Madau (2012) at redshift $z=0$.

In the top panel of Figure 2, we show the gas temperature and the mean molecular weight $\mu$ as a function of a single timestep $\Delta t$ for a starting temperature of $10\,\rm{K}$. The gas temperature is expected to rise up to \sim$$10^{4}\,\rm{K} due to the heating effect of the UV background, which completes at $\sim$1 Myr. The predicted temperature of the metal-free gas case is very similar to $0.1Z_{\odot}$ gas.

The middle panel shows the gas temperature and $\mu$ for a starting temperature of $10^{6}$ K, which is the same as that in Smith et al. (2017). The gas temperature and $\mu$ for $0.1Z_{\odot}$ gas agree perfectly with Figure 3 in Smith et al. (2017). In the case of net cooling, metal-free gas cools more slowly than $0.1Z_{\odot}$, which is also subject to metal cooling. Note that cooling from $10^{5}$ to $10^{4}$ K is a very abrupt phase, much shorter than what it takes from $10^{3}$ to $10^{4}$ 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 $\Delta t\approx 15\,\rm{Myr}$, 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 zoomed-in view of the gas temperature and $\mu$ 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 $\mu$ accurately reflect the features in the cooling curves between $10^{4}$ and $10^{5}$ K (a local minimum between two peaks) within a short duration of $\sim$0.2 Myr.

3.2 Cosmological hydrodynamic simulations with different time integration schemes for gas cooling

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 $\Lambda$CDM using the music code (Hahn & Abel, 2011). A uniform box of $30$ comoving Mpc on each side contains a total number of $256^{3}$ gas particles and $256^{3}$ dark matter particles. Hydrodynamics is evolved using the moving-finite-mass (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 exact-integration 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 Faucher-Giguè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 $z=0.6$ simulated with the cooling method in Gadget, Grackle, and exact integration scheme in the entire simulation box. Gas temperature is color-coded 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 zoomed-in 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 low-density 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 density-temperature 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 accurate444Since 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 $z=2$ and $z=0.6$, the equilibrium between heating and cooling drops from \sim$$10^{5}\,\rm{K} for gas at density $n_{\rm H}=10^{-7}\,\rm{cm^{-3}}$ to \sim$$10^{4}\,\rm{K} for $n_{\rm H}=10^{-1}\,\rm{cm^{-3}}$.

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 $10^{-2}\leq n_{\rm H}\leq 10^{-1}\,\rm{cm^{-3}}$ is slightly higher in the exact time integration than the implicit integration and Grackle, which contains many gas particles with a temperature below $10^{4}$ 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 $z=1.1$ and $z=0.6$ in $10^{-2}\leq n_{\rm H}\leq 10^{-1}\,\rm{cm^{-3}}$. In particular, the distribution of the gas temperature is skewed below $10^{4}$ K for gas at densities around $10^{-1}\,\rm{cm^{-3}}$. Interestingly, the gas temperature follows the equilibrium relation very well at $z=2.2$. For the simulation using Grackle cooling, the gas temperature is consistently lower than the equilibrium line in the density range $10^{-2}\leq n_{\rm H}\leq 10^{-1}\,\rm{cm^{-3}}$.

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 $10^{5}$ to $\sim 10^{4}$ K than heat it from $10^{3}$ to $\sim 10^{4}$ 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 $M_{\rm DM}$–$M_{\rm stellar}$ 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 over-cooling 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 $z=0.6$ 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 $M_{\rm stellar}^{\rm exact}$ as a baseline for comparison. The differences in stellar mass is shown in the lower part of Figure 5 in terms of $\Delta\log(M_{\rm stellar})=\log(M_{\rm stellar})-\log(M_{\rm stellar}^{\rm exact})$. The stellar mass for a given dark matter halo mass in the three simulations agrees to within $0.1$ dex. Nevertheless, the differences in stellar mass between the exact approach and the other integrators is systematic for dark matter haloes more massive than $10^{11}M_{\rm\odot}$, with the exact technique yielding the lowest masses. Since the peak of star formation efficiencies occurs at \sim$$10^{12}M_{\rm\odot} 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 well-defined 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 $10^{5}$ yr in Wiersma et al. (2009) indicates that there could be certain gas phases not well-resolved 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 $\Delta t$ 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 $\Delta t$ for each particle at $z=0.6$ with the exact time integration. We use a simple mapping from the ratios to colors of gas particles as shown in the color-bar. A very long cooling time with respect to $\Delta t$ 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 shock-heated region, the low density gas with temperature above $10^{6}$ 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 phase-space 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 $\log(T)\sim 5$ and $n_{\rm H}\sim 10^{-7}\,\rm{cm^{-3}}$. For gas particles in the filaments ($\log(T)\sim 4$ and $10^{-2}\leq n_{\rm H}\leq 10^{-1}\,\rm{cm^{-3}}$), the gas temperature lies around the equilibrium in the exact integration scheme. On the other hand, $t_{\rm cool}/\Delta t$ in this phase does not exceed $10^{5}$ as the green line terminates around $n_{\rm H}=10^{-4}\,\rm{cm^{-3}}$. 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 $t_{\rm cool}/\Delta t\sim 0.01$. 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. Ultra-high 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 self-shielding (Rahmati et al., 2013), or a local radiation field (Vogelsberger et al., 2013; Gnedin & Hollon, 2012) or molecular hydrogen ($\rm{H_{2}}$) cooling, which requires one to generate cooling tables using photo-ionization codes such as CLOUDY. On the other hand, this method is difficult to apply to non-equilibrium cooling.

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 wall-time statistics of the three integration schemes used in our simulations. The wall-time 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 look-ups and interpolations, which we have not optimized, during its root finding step. In theory, the cost of the exact integration scales as $N_{\rm gas}N_{\rm grid}$, the product of the number of gas particles $N_{\rm gas}$ and the number of temperature grid points $N_{\rm grid}$.

This time integration scheme might also help to resolve gas temperature for studies of cooling emission (Haiman et al., 2000; Fardal et al., 2001; Faucher-Giguère et al., 2010b). Because gas at temperatures around \sim$$10^{4} K in the current simulations is not well-resolved 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-$\alpha$ emissivity changes exponentially around \sim$$10^{4} K (Faucher-Giguère et al., 2010a).

While the differences in the stellar mass produced by the three time integration schemes for gas cooling are relatively small ($<0.1$ 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 $0.1\,\rm{cm^{-3}}$) 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 well-resolved except for gas more dense than $n_{\rm H}\sim 0.01\ {\rm cm^{-3}}$. 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 redshift-dependent 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

1. (1)

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. 2. (2)

It is easy to implement and efficient to execute. 3. (3)

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

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 AST-0965694, AST-1009867, and AST-1412719. 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.