The Overcooling Problem and the Sub-Grid Recipe: Where Galaxy Simulations Stop Being First-Principles — Epoche C2
This note sets out a problem that has been the organising difficulty of numerical galaxy formation for four decades, the two-part solution the field has settled on, and the reason that solution does not make the simulations first-principles calculations. Units follow the astrophysical literature throughout: masses in $M_\odot$, distances in pc and kpc, speeds in km s$^{-1}$, energies in erg, number densities in cm$^{-3}$. 1. The problem: cooling is too efficient In the White-Rees picture, gas falls into a dark matter halo, shocks to the virial temperature, radiates, and condenses. The question is how fast. The cooling time is the thermal energy per unit volume divided by the radiated power per unit volume, $$t_{\mathrm{cool}} \simeq \frac{3 n_H k_B T}{n_H^2 \Lambda(T)} = \frac{3 k_B T}{n_H \Lambda(T)},$$ where $n_H$ is the hydrogen number density, $T$ the gas temperature, $k_B = 1.381\times10^{-16}$ erg K$^{-1}$, and $\Lambda(T)$ the cooling function — the power radiated per unit volume divided by the square of the density, so that it depends on temperature and composition but not on how much gas there is. The factor 3 rather than 3/2 accounts for the free electrons in ionised gas. Evaluate this for the hot halo of a Milky Way-like galaxy. Take $T = 10^{6}$ K, $n_H = 10^{-3}$ cm$^{-3}$, and $\Lambda = 10^{-23}$ erg cm$^{3}$ s$^{-1}$, a representative value near $10^6$ K. Then $3k_BT = 4.1\times10^{-10}$ erg and $n_H\Lambda = 10^{-26}$ erg s$^{-1}$, giving $t_{\mathrm{cool}} = 4.1\times10^{16}$ s $= 1.3$ Gyr. Compare the halo dynamical time: with $R_{\mathrm{vir}} = 200$ kpc $= 6.2\times10^{23}$ cm and circular speed $v_c = 200$ km s$^{-1} = 2\times10^{7}$ cm s$^{-1}$, $t_{\mathrm{dyn}} = R_{\mathrm{vir}}/v_c = 3.1\times10^{16}$ s $= 1.0$ Gyr. The two are equal to within 30%. Gas cools as fast as it falls. A simulation containing only gravity, hydrodynamics and radiative cooling therefore converts most available baryons into stars. Reality does not. The cosmic baryon fraction is $f_b = \Omega_b h^2/\Omega_m h^2 = 0.0224/0.143 = 0.156$. For a halo of $1.5\times10^{12}\,M_\odot$ the baryon budget is $2.3\times10^{11}\,M_\odot$, against a Milky Way stellar mass of $5\times10^{10}\,M_\odot$ — an efficiency of $0.21$. And this is the maximum. At $M_h = 10^{10}\,M_\odot$ a typical dwarf has $M_\star \sim 10^{7}\,M_\odot$, an efficiency of $10^{7}/(0.156\times10^{10}) = 0.0064$, or 0.6%. The efficiency falls on both sides of a peak near the Milky Way's halo mass, by a factor of thirty within two decades of halo mass. 2. The solution, part one: supernovae, and why they stop working A Chabrier initial mass function yields roughly one core-collapse supernova per $100\,M_\odot$ of stars formed, each releasing $E_{\mathrm{SN}} \approx 10^{51}$ erg. The specific energy budget is therefore $10^{49}$ erg per $M_\odot$ of stars. Ask how much gas that can eject. For an energy-driven wind, $\tfrac{1}{2}\eta M_\star v_{\mathrm{esc}}^2 = f \times 10^{49}\,\mathrm{erg}\,M_\odot^{-1} \times M_\star$, where $\eta$ is the mass loading (gas ejected per unit stellar mass formed) and $f$ the fraction of supernova energy that survives radiative losses. Taking $v_{\mathrm{esc}} \approx 2v_c$: Dwarf, $v_c = 50$ km s$^{-1}$: $v_{\mathrm{esc}} = 10^{7}$ cm s$^{-1}$, so $v_{\mathrm{esc}}^2 = 10^{14}$ erg g$^{-1} = 2.0\times10^{47}$ erg per $M_\odot$. With $f = 0.1$, $\eta = 2\times0.1\times10^{49}/2.0\times10^{47} = 10$. Milky Way, $v_c = 200$ km s$^{-1}$: $v_{\mathrm{esc}}^2 = 3.2\times10^{48}$ erg per $M_\odot$, giving $\eta = 0.63$. The scaling is $\eta \propto v_c^{-2}$, and the ratio of the two cases is $(200/50)^2 = 16$, which is what the numbers show. Supernovae expel ten times their own stellar mass from a dwarf and cannot lift their own weight out of a massive galaxy. This is the first half of the answer, and it explains the low-mass side of the peak only. 3. The solution, part two: the black hole For the high-mass side the available reservoir is accretion. A black hole of mass $M_{\mathrm{BH}}$ radiates $E_{\mathrm{AGN}} = \varepsilon M_{\mathrm{BH}}c^2$ with $\varepsilon \approx 0.1$. Observed black holes satisfy $M_{\mathrm{BH}} \approx 1.4\times10^{-3}M_{\mathrm{bulge}}$. The binding energy of the bulge is of order $M_{\mathrm{bulge}}\sigma^2$, with $\sigma$ the stellar velocity dispersion. The ratio is $$\frac{E_{\mathrm{AGN}}}{E_{\mathrm{bind}}} = \frac{\varepsilon\,(1.4\times10^{-3})\,c^2}{\sigma^2} = \frac{1.4\times10^{-4}\times8.99\times10^{20}}{4\times10^{14}} \approx 3\times10^{2},$$ using $\sigma = 200$ km s$^{-1}$, so $\sigma^2 = 4\times10^{14}$ erg g$^{-1}$. This is the Silk-Rees argument. Even a coupling efficiency of one per cent leaves a factor of three in hand, exactly where supernovae have run out of energy. The two mechanisms are weakest in the same place, and that place is where the stellar-to-halo mass relation peaks. 4. Why this is a recipe and not a calculation The physics above happens on scales the simulations do not resolve. A supernova remnant remains energy-conserving until it has swept up of order $10^{3}\,M_\odot$; at $n_H = 1$ cm$^{-3}$ that is a radius of about 30 pc, since $\tfrac{4}{3}\pi(9.3\times10^{19}\,\mathrm{cm})^3 \times 1\,\mathrm{cm}^{-3}\times 1.67\times10^{-24}$ g $= 2.8\times10^{3}\,M_\odot$. After that the remnant radiates most of its energy and only momentum survives. Large cosmological volumes reach a baryonic mass resolution of order $10^{6}\,M_\odot$ and a gravitational softening of about 0.7 kpc. The shortfalls are $10^{6}/10^{3} = 10^{3}$ in mass, three orders of magnitude, and $700/30 = 23$ in length, between one and two. The consequence is quantitative. Deposit $10^{51}$ erg as heat into a single gas element of $10^{6}\,M_\odot$, which contains $1.99\times10^{39}/1.67\times10^{-24} = 1.2\times10^{63}$ particles. The temperature rise is $\Delta T = \tfrac{2}{3}E/(Nk_B) = \tfrac{2}{3}\times10^{51}/(1.2\times10^{63}\times1.38\times10^{-16}) = 4\times10^{3}$ K