Cooling flow solutions for the circumgalactic medium

TL;DR

This paper develops analytic and simulation-based models of cooling flows in galaxy halos, showing they can explain certain observations and serve as benchmarks for feedback processes in galaxy formation.

Contribution

It introduces a family of steady-state cooling flow solutions for galaxy halos, linking inflow rates to halo properties and matching several observational features.

Findings

Cooling flows can predict observed absorption and emission features.

Steady-state solutions are determined by a single parameter, such as inflow rate.

Some galaxy clusters' entropy profiles align with cooling flow models.

Abstract

In several models of galaxy formation feedback occurs in cycles or mainly at high redshift. At times and in regions where feedback heating is ineffective, hot gas in the galaxy halo is expected to form a cooling flow, where the gas advects inward on a cooling timescale. Cooling flow solutions can thus be used as a benchmark for observations and simulations to constrain the timing and extent of feedback heating. Using analytic calculations and idealized 3D hydrodynamic simulations, we show that for a given halo mass and cooling function, steady-state cooling flows form a single-parameter family of solutions, while initially hydrostatic gaseous halos converge on one of these solutions within a cooling time. The solution is thus fully determined once either the mass inflow rate $\dot{M}$ or the total halo gas mass are known. In the Milky Way (MW) halo, a cooling flow with $\dot{M}$ equal to the star formation rate predicts a ratio of the cooling time to the free-fall time of ~10, similar to some feedback-regulated models. This solution also correctly predicts observed OVII and OVIII absorption columns, and the gas density profile implied by OVII and OVIII emission. These results suggest ongoing heating by feedback may be negligible in the inner MW halo. Extending similar solutions out to the cooling radius however underpredicts observed OVI columns around the MW and around other low-redshift star-forming galaxies. This can be reconciled with the successes of the cooling flow model with either a mechanism which preferentially heats the OVI-bearing outer halo, or alternatively if OVI traces cool photoionized gas beyond the accretion shock. We also demonstrate that the entropy profiles of some of the most relaxed clusters are reasonably well described by a cooling flow solution.

Equations

$$\dot{M}=4\pi r^2\rho v$$

$$\frac{1}{2}\frac{\mathrm{d}{v}^2}{\mathrm{d}r}=-\frac{1}{\rho }\frac{\mathrm{d}P}{\mathrm{d}r}-\frac{v_{\mathrm{c}}^2}{r}$$

$$v\frac{\mathrm{d}\ln K}{\mathrm{d}r}=-\frac{1}{t_{\mathrm{cool}}}$$

$$v_{\mathrm{c}}\equiv \sqrt{\frac{G{M}_{\mathrm{grav}}(<r)}{r}}$$

$$t_{\mathrm{cool}}=\frac{P}{(\gamma -1)n_{\mathrm{H}}^2\Lambda }.$$

$$t_{\mathrm{ff}}=\frac{\sqrt{2}r}{v_{\mathrm{c}}}$$

$$\mathcal{B}\equiv \frac{v^2}{2}+\frac{c_{\mathrm{s}}^2}{\gamma -1}+\Phi$$

$$R_{\mathrm{vir}}=260{\left(\frac{M_{\mathrm{halo}}}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{1/3}\mathrm{kpc}$$

$$T_{\mathrm{vir}}=\frac{\mu m_{\mathrm{p}}{v}_{\mathrm{c}}(R_{\mathrm{vir}}{)}^2}{2k}=6\cdot 10^5{\left(\frac{M_{\mathrm{halo}}}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{2/3}\mathrm{K}.$$

$$M_{\ast }(<R)=M_{\ast }\frac{r}{r+R_{1/2}}$$

$$R_{1/2}=0.015R_{200c}=3.0{\left(\frac{M_{\mathrm{halo}}}{10^{12}{\mathrm{M}}_{\odot }}\right)}^{1/3}\mathrm{kpc}.$$

$$\frac{\mathrm{d}\ln \rho }{\mathrm{d}\ln r}+\frac{\mathrm{d}\ln v}{\mathrm{d}\ln r}=-2,$$

$$\frac{\mathrm{d}\ln T}{\mathrm{d}\ln r}-(\gamma -1)\frac{\mathrm{d}\ln \rho }{\mathrm{d}\ln r}=\frac{t_{\mathrm{flow}}}{t_{\mathrm{cool}}},$$

$$t_{\mathrm{flow}}\equiv r/|v|.$$

$$c_{\mathrm{s}}=\sqrt{\frac{\gamma P}{\rho }}$$

$${\mathcal{M}}^2\frac{\mathrm{d}\ln v}{\mathrm{d}\ln r}=-\frac{1}{\gamma }\frac{\mathrm{d}\ln P}{\mathrm{d}\ln r}-\frac{v_{\mathrm{c}}^2}{c_{\mathrm{s}}^2}$$

$$\frac{\mathrm{d}\ln v}{\mathrm{d}\ln r}\left({\mathcal{M}}^2-1\right)=2-\frac{v_{\mathrm{c}}^2}{c_{\mathrm{s}}^2}-\frac{t_{\mathrm{flow}}}{\gamma t_{\mathrm{cool}}}.$$

$$\max (R_{1/2},R_{\mathrm{feedback}})\lesssim r\lesssim \min (R_{\mathrm{cool}},R_{\mathrm{shock}})$$

$$v_{\mathrm{c}}(r)=v_{\mathrm{c}}(R_{\mathrm{vir}}){\left(\frac{r}{R_{\mathrm{vir}}}\right)}^m$$

$$T\propto r^{2m},n_{\mathrm{H}}\propto r^{-\frac{3}{2}+m},\mathcal{M}\propto r^{-\frac{1}{2}-2m}$$

$$\frac{v_{\mathrm{c}}^2}{c_{\mathrm{s}}^2}=\frac{9}{10}(1-2m)\equiv A,$$

$$\frac{t_{\mathrm{flow}}}{t_{\mathrm{cool}}}=1+\frac{4}{3}m\equiv B.$$

$$T_{\mathrm{c}}(r)\equiv \frac{\mu m_{\mathrm{p}}{v}_{\mathrm{c}}^2(r)}{\gamma k_{\mathrm{B}}}$$

$$T(r)=\frac{1}{A}{T}_{\mathrm{c}}(r)=\frac{6}{5A}{\left(\frac{r}{R_{\mathrm{vir}}}\right)}^{2m}{T}_{\mathrm{vir}}.$$

$$K\propto r^B\propto r^{1+4m/3}.$$

$$\mathcal{M}=\frac{XA}{m_{\mathrm{p}}{v}_{\mathrm{c}}^2}\sqrt{\frac{5\dot{M}\Lambda }{18\pi Br}}$$

$$n_{\mathrm{H}}=\sqrt{\frac{9B\dot{M}}{40\pi A{r}^3\Lambda }}{v}_{\mathrm{c}}.$$

$$\frac{t_{\mathrm{cool}}}{t_{\mathrm{ff}}}=\frac{t_{\mathrm{flow}}/B}{\sqrt{2}r/v_{\mathrm{c}}}=\frac{\sqrt{A}}{\sqrt{2}B}{\mathcal{M}}^{-1},$$

$$n_{\mathrm{H}}=1.6\cdot 10^{-5}{M}_{12}^{0.36}{\dot{M}}_1^{1/2}{\Lambda }_{-22}^{-1/2}{r}_{100}^{-1.6}{\mathrm{cm}}^{-3},$$

$$t_{\mathrm{cool}}=7.2M_{12}^{0.36}{\dot{M}}_1^{-1/2}{\Lambda }_{-22}^{-1/2}{r}_{100}^{1.4}\mathrm{Gyr},$$

$$R_{\mathrm{cool}}=130M_{12}^{-0.26}{\dot{M}}_1^{0.36}{\Lambda }_{-22}^{0.36}\mathrm{kpc}.$$

Full text

††pubyear: 2019††pagerange: Cooling flow solutions for the circumgalactic medium–C

1 Introduction

Classic cooling flow solutions were derived in the 1980’s following X-ray observations of cluster centers, which revealed gas cooling times shorter than the Hubble time. In these spherical solutions the loss of entropy via radiation drives an inflow towards the center of the potential well, at a rate where heating by compression roughly balances radiative energy losses so that the gas temperature remains near the virial temperature (see reviews in Fabian et al. 1984; Sarazin 1986). The inflowing gas is expected to be predominantly single phase, since the flow is advected on the same timescale as the timescale on which thermal instabilities grow (Mathews & Bregman 1978; Malagoli et al. 1987; Balbus & Soker 1989; Li & Bryan 2012). Cooling flow solutions have however fallen out of favor in the context of clusters, since the predicted mass flow rate based on the observed X-ray emission exceeds the observed star formation rate in the central galaxy by a factor of $10-100$, and the models predict a central spike in emission which is not observed (see reviews by McNamara & Nulsen 2007; Fabian 2012; McDonald et al. 2018).

Testing similar solutions in halos less massive than clusters has been hampered by the lower emission measure of the virial temperature gas, which drops sharply with decreasing halo mass. Studies have hence been mostly limited to comparing observed star formation rates (SFRs) with estimates of cooling mass flow rates ${\dot{M}}$ based on some assumption on the mass of the cooling circumgalactic gas, say that it is equal to the halo baryon budget (e.g. White & Frenk 1991). In these studies ${\dot{M}}$ typically overpredicts the observed SFRs. However, since in a cooling flow ${\dot{M}}\propto M_{\rm gas}^{2}$, where $M_{\rm gas}$ is the cooling gas mass, if $M_{\rm gas}$ is overestimated then ${\dot{M}}$ will be overpredicted. The value of $M_{\rm gas}$ could potentially be lower than the halo baryon budget especially at low redshift, due to the effects of strong winds at earlier epochs (e.g., Muratov et al. 2015; Hafen et al. 2018). Moreover, most of the cooling gas may be ejected from the galaxy as outflows rather than form stars, in which case the SFRs underestimate ${\dot{M}}$.

In recent years deep X-ray observations with XMM-Newton and Chandra have detected absorption and emission from O vii and O viii lines, which most likely originate in the hot gaseous halo of the Milky Way (Henley & Shelton 2010, 2012; Gupta et al. 2012; Miller & Bregman 2013; Fang et al. 2015; Bregman et al. 2018). The installation of the Cosmic Origins Spectrograph (COS) onboard HST also facilitated surveys of high-ionization ions such as O vi and Ne viii in sightlines through halos of external galaxies (Chen & Mulchaey 2009; Prochaska et al. 2011; Tumlinson et al. 2011; Johnson et al. 2015, 2017; Keeney et al. 2018; Chen et al. 2018; Burchett et al. 2019). Unless the absorbing gas pressure is a factor of $\sim 30$ less than expected in a virial-temperature gaseous halo, these UV absorption features must also trace gas at or near the virial temperature (McQuinn & Werk 2018; Stern et al. 2018; Burchett et al. 2019). Thus, using these observational constraints, we can test cooling flow solutions in halos less massive than clusters without any assumption on the cooling gas mass or the fraction of it which turns into stars.

Cooling flow solutions are also useful as a benchmark for cosmological and idealized hydrodynamic simulations, in order to disentangle the effects of radiative cooling on gaseous halos from the effects of other physical processes, such as galaxy feedback. Previous studies have typically compared gaseous halos in simulations to hydrostatic solutions (e.g. McCarthy et al. 2010; Fielding et al. 2017; Oppenheimer 2018). Since in hydrostatic solutions cooling is neglected, it is not straightforward to discern whether differences between the simulation and the idealized solution are due to cooling or due to other physical effects. Cooling flow solutions bypass this limitation by providing the expected physical properties of gaseous halos if only cooling is present.

Cooling flows are also an integral phase of time-dependent models of feedback such as feedback limit cycles, in which a cooling flow initially develops, fuels feedback which in turn suppresses the inflow, after which the supply of gas needed to maintain feedback stops and a cooling flow redevelops, and so forth (e.g. Ciotti & Ostriker 2001; Pizzolato & Soker 2005; Li & Bryan 2014a, b; Meece et al. 2015; Prasad et al. 2015, 2017; Soker 2016; Yang & Reynolds 2016; Martizzi et al. 2019). Understanding cooling flows and their observable consequences is hence useful to test this class of models, and to constrain the duty-cycle of their cooling flow phase.

The goal of this paper is to systematically adapt the cooling flow solutions originally derived for cluster-scale halos, to a wide range of halo masses, down to the scale of dwarf galaxy halos. We include in the solutions updated constraints on the distribution of matter and properties of radiative cooling, which were not available when the original cooling flow solutions where derived. We then proceed to test these cooling flow solutions against the recent X-ray and UV constraints mentioned above. We focus on redshift $z\sim 0$ halos where observations are most constraining, and defer an analysis of halos at higher redshifts to future work. In two companion papers we discuss the implications of cooling flow solutions for the transition between ‘cold mode’ accretion (=‘cold flows’) and ‘hot mode’ accretion (=‘cooling flows’), and compare cooling flow solutions with gaseous halos in the FIRE-2 cosmological simulations (Hopkins et al. 2018). These papers are referred to below as Paper II and Paper III, respectively.

This paper is organized as follows. In section 2 we find cooling flow solutions to the steady-state flow equations, first by deriving a self-similar solution akin to Fabian et al. (1984), and then by direct integration of the flow equations. In section 3 we test the validity of these steady-state solutions using idealized 3D hydrodynamic simulations. Section 4 compares the predictions of cooling flow solutions with available observational constraints. We discuss and compare our results to previous work in section 5 and summarize in section 6.

A flat $\Lambda$CDM cosmology with $H_{0}=68\,\rm km\ s^{-1}\,{\rm Mpc}^{-1}$, $\Omega_{\rm M}=0.31$, and a cosmic baryon fraction $f_{\rm b}=0.158$ is assumed throughout (Planck Collaboration et al. 2016).

2 Steady-state equations

In this section we solve the spherical steady-state equations for cooling flows. We first present definitions and equations (section 2.1), and discuss the range of radii where our simplifying assumptions could be applicable (section 2.2). We then find self-similar solutions (section 2.3) and full numerical solutions to the flow equations (section 2.4). In the last subsections 2.5–2.6 we discuss the properties of the solutions for specific halo masses.

2.1 Equations and definitions

The steady-state equations for mass, momentum, and entropy conservation of a spherically-symmetric ideal fluid, with no angular momentum, magnetic field, viscosity, or thermal conduction, are

[TABLE]

In these equations $r$, $v$, $\rho$, $P$, and $\ln K$ are respectively the radius, radial velocity (negative for an inflow), gas density, gas pressure, and gas entropy ($K\propto P/\rho^{\gamma}$, where $\gamma$ is the adiabatic index). Also,

[TABLE]

is the circular velocity where $M_{\rm grav}(<r)$ is the gravitating mass within $r$, and $t_{\rm cool}$ is the cooling time, defined as the ratio of the energy per unit volume $(\gamma-1)^{-1}P$ to the radiated energy per unit volume $n_{\rm H}^{2}\Lambda$:

[TABLE]

where $n_{\rm H}=X\rho/m_{\rm p}$ is the hydrogen density ($X$ is the hydrogen mass fraction) and $\Lambda$ is the cooling function. We also define the free-fall time

[TABLE]

and the Bernoulli parameter

[TABLE]

where $\Phi=-\int(v_{\rm c}^{2}/r){\rm d}r$ is the gravitational potential.

The left panel of Figure 1 plots the values of $v_{\rm c}(r)$ used in this paper for redshift $z=0$ halos with different masses. The assumed radial mass distributions are composed of three components: an NFW halo, a central galaxy, and an outer component. We ignore the contribution of gas mass to the gravitational mass, validating this assumption a postriori with our derived gas models. For the NFW halo, we calculate the concentration parameter using the relation in Dutton & Macciò (2014). The implied virial radius for an overdensity defined as in Bryan & Norman (1998) is

[TABLE]

(we use capital $R$ for radial quantities and small $r$ for the coordinate), while the virial temperature is

[TABLE]

The central galaxy is assumed to have a stellar mass $M_{*}$ estimated from $M_{\rm halo}$ based on the Behroozi et al. (2018) relation for central galaxies. For simplicity, the stellar mass distribution is assumed to have a radial profile:

[TABLE]

with the stellar half-mass radius $R_{1/2}$ taken from Kravtsov (2013):

[TABLE]

Here, $R_{\rm 200c}$ is the radius enclosing an overdensity of 200 relative to the critical density. The values of $R_{1/2}$ are marked in the left panel of Fig. 1. Individual galaxies have a scatter of $0.2\,{\rm dex}$ around this value (Kravtsov 2013). The value of $R_{1/2}$ is also an estimate of the radius where angular momentum may become important (see below). For the outer halo component we use the formulation deduced by Diemer & Kravtsov (2014, hereafter DK14), with their median values of $s_{\rm e}=1.5$ and $b_{\rm e}=1$.

The right panel of Fig. 1 plots $\Lambda(T,n_{\rm H},Z)$ in $z=0$ halos ($T$ and $Z$ are the gas temperature and metallicity, respectively), based on the calculations of Wiersma et al. (2009), which assume optically thin ionization equilibrium conditions and the UV background of Haardt & Madau (2012, hereafter HM12). When calculating $\Lambda$ we use a density which is half the assumed density, since the HM12 background likely underestimates the background by a factor of two (Faucher-Giguère 2019).

We show below that it is convenient to use the logarithmic derivatives of $\rho$, $v$ and $T$, so we cast the flow equations (1)–(3) in logarithmic form. The mass equation (1) is equivalent to

[TABLE]

while the entropy equation (3) can be cast as

[TABLE]

where we used the definition of the flow time

[TABLE]

Also, defining the adiabatic sound speed

[TABLE]

and multiplying the momentum eqn. (2) by $r/c_{\rm s}^{2}$ yields:

[TABLE]

where $\mathcal{M}\equiv|v|/c_{\rm s}$ is the Mach number. Using eqns. (12) and (13) to cancel ${\rm d}\ln P={\rm d}\ln T+{\rm d}\ln\rho$, we get after some rearranging

[TABLE]

For $t_{\rm cool}\rightarrow\infty$ equation (17) reduces to the standard equation used in analyzing adiabatic Bondi flows.

2.2 Radii where solutions are applicable

Before finding solutions to the steady-state equations (1)–(3), we note that the maximum range of radii where we expect them to apply is:

[TABLE]

where $R_{1/2}$ is used here as an approximation of the radius where the halo gas could be supported against gravity by angular momentum, $R_{\rm feedback}$ is a putative maximum radius where feedback by the galaxy heats, or otherwise alters the physical properties of the circumgalactic gas, $R_{\rm cool}$ is the usual ‘cooling radius’, where the cooling time equals the age of the system or the time since the last heating event, and $R_{\rm shock}$ is the radius of the accretion shock. We now explain each of these terms.

The $R_{1/2}$ limit appears on the left side of eqn. (18) since the flow equations neglect angular momentum, while centrifugal forces will be significant on the galaxy scale if the specific angular momentum of the gas is similar to the average of the dark matter halo (e.g. Kravtsov 2013). The feedback radius $R_{\rm feedback}$ is hard to estimate a priori, and one of the main goals of this study is to provide a benchmark solution in which $R_{\rm feedback}=0$, so deviations from this solution in observations and simulations could be used to constrain $R_{\rm feedback}$. The $R_{\rm cool}$ limit on the right-hand side of eqn. (18) is because beyond this radius a steady-state cooling flow does not have time to develop, and a time-dependent solution such as Bertschinger (1989) must be found. The value of $R_{\rm cool}$ is derived below as a function of halo gas properties. Last, a discontinuity is expected at $R_{\rm shock}$, beyond which the gas is free-falling and close to thermal equilibrium with the UV background, so a cooling flow solution does not apply.

2.3 Subsonic self-similar solution

Following Fabian et al. (1984), we derive a self-similar solution to the flow equations in the subsonic limit $(\mathcal{M}^{2}\rightarrow 0)$, by approximating $v_{\rm c}$ as a power-law:

[TABLE]

where $m$ is constant. For an isothermal potential $m=0$, while around a point mass $m=-0.5$. The self-similar solution can then be found by requiring that all logarithmic derivatives of the gas properties are constant, in which case the ratios $t_{\rm cool}/t_{\rm flow}$ and $v_{\rm c}^{2}/c_{\rm s}^{2}$ are also constant (see eqns. 13 and 17). For constant $\Lambda$, these two conditions yield (see above definitions of $t_{\rm cool}$, $t_{\rm flow}$, and $c_{\rm s}$):

[TABLE]

where here and henceforth we use $\gamma=5/3$. More general relations can similarly be found for any power-law dependence of $\Lambda$ on $T$, $n_{\rm H}$, or $r$, though we avoid this complication since it only mildly increases the accuracy of our analytic estimates while significantly increasing the complexity of the analytic expressions. In the numerical integration in the next section we use the full forms of $v_{\rm c}$ and $\Lambda$ plotted in Fig. 1.

Plugging the relations in eqn. (20) into equations (13) and (17) yields

[TABLE]

and

[TABLE]

Eqn. (21) implies that the gas temperature in cooling flows is approximately equal to the ‘circular temperature’, defined as

[TABLE]

which is noted in the left panel of Fig. 1. We hence get

[TABLE]

where in the second equality we used the definition of $T_{\rm vir}$ (eqn. 9). Also, eqn. (22) together with eqn. (13) imply that the entropy profile in cooling flows scales as

[TABLE]

This proportionality was previously derived by McCarthy et al. (2005), and can also be derived from eqn. (12) in Voit (2011), who found $K\propto rT^{2/3}$ by assuming $t_{\rm cool}\approx t_{\rm flow}$ and constant ${\dot{M}}$. Adding the self-similar requirement that $T\propto r^{2m}$ (eqn. 20) then yields eqn. (25).

To complete the solution, we need to derive also the density and Mach number, which depend on the free parameter ${\dot{M}}$. Extracting $\mathcal{M}$ from ${\dot{M}}=4\pi r^{2}\rho v$ (eqn. 1) and $v=r/Bt_{\rm cool}$ (eqn. 22) we get

[TABLE]

where we used eqn. (5) for $t_{\rm cool}$, eqn. (21) for $c_{\rm s}^{2}$, and we remind the reader that $X$ is the hydrogen mass fraction. Similarly, extracting $n_{\rm H}$ gives

[TABLE]

From equations (24)–(27) and the numerical values of $A$ and $B$ given in eqns. (21)–(22), all the properties of the solution can be found. These equations show that for a given $M_{\rm halo}$ and $\Lambda$, cooling flow solutions have a single free parameter. This free-parameter can be chosen according to convenience, and below we use either ${\dot{M}}$, the total gas mass $M_{\rm gas}$, or the sonic radius $R_{\rm sonic}$ where $\mathcal{M}=1$.

An important property of cooling flows is the relation between $t_{\rm cool}/t_{\rm ff}$ and $\mathcal{M}$. From eqns. (21)–(22) and the definition of $t_{\rm ff}$ in eqn. (6) we get

[TABLE]

i.e. in cooling flows $t_{\rm cool}/t_{\rm ff}$ is equal to $\mathcal{M}^{-1}$ up to a factor of order unity. At $R_{\rm sonic}$ we expect $t_{\rm cool}\approx t_{\rm ff}$.

Below we use these self-similar solutions to derive analytic estimates of the physical properties of cooling flows. We first though compare these solutions to more accurate solutions derived by direct integration.

2.4 Integration of flow equations

Direct integration of the flow equations allows us to account for the effects of a finite $\mathcal{M}$ on the solution, and also incorporates the full forms of $v_{\rm c}$ and $\Lambda$ shown in Fig. 1, rather than the approximations used in the previous section. Similar integrations have previously been done by Cox & Smith (1976) and Mathews & Bregman (1978). We start the integration from a sonic point at some radius $R_{\rm sonic}$ and integrate outward, allowing $R_{\rm sonic}$ to be outside the range of radii of interest (i.e. not all solutions are transonic at halo radii). At $R_{\rm sonic}$ the right side of eqn. (17) must vanish or the velocity derivative will be infinite, hence the three boundary values $\rho(R_{\rm sonic})$, $T(R_{\rm sonic})$ and $v(R_{\rm sonic})$ must satisfy the two conditions $\mathcal{M}=1$ and $2-v_{\rm c}^{2}/c_{\rm s}^{2}-t_{\rm flow}/\gamma t_{\rm cool}=0$. A third condition on the boundary values is deduced by requiring that the solution is marginally bound, i.e. that $\mathcal{B}\rightarrow 0^{-}$ as $r\rightarrow\infty$. This choice of the outer boundary condition is motivated by the resemblance to the self-similar solutions discussed in the previous section, which are also marginally-bound. This boundary condition is enforced on the solution via a shooting method described in Appendix A, and yields a tight constraint on $T(R_{\rm sonic})$. There is thus only a single marginally-bound solution for a given $R_{\rm sonic}$, and the entire process yields a one-parameter family of solutions which differ in their value of $R_{\rm sonic}$.

The details of the integrating through the sonic point are also given in Appendix A. We note that other choices of the outer boundary condition, such as an accretion shock or a specific temperature at the cooling radius, change the solution relative to the marginally-bound solution only near the outer boundary, and hence do not affect our conclusions. This insensitivity of the solutions to the exact choice of the boundary condition is demonstrated in Appendix B.

In Figure 2 we plot an example transonic marginally-bound solution, for $M_{\rm halo}=10^{12}\,{\rm M_{\odot}}$, $Z=0.3\,{\rm Z_{\odot}}$, and an assumed $R_{\rm sonic}=1\,{\rm kpc}$. The blue lines in the different panels show different properties of the solution. Going inward the temperature roughly follows $T_{\rm c}$ (top-left panel), while the entropy drops (third panel on the left). This behavior demonstrates that energy loss to radiation is compensated by adiabatic compression, and hence the effect of radiative cooling is most apparent as a drop in entropy rather than a drop in temperature. Gray lines in the panels plot a self-similar solution with the same ${\dot{M}}=4.5\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ as the integrated solution, using the approximations $\Lambda_{-22}=0.6$ and $v_{\rm c}=140r_{100}^{m}\,\rm km\ s^{-1}$ with $m=-0.1$, where $r=100\,r_{100}\,{\rm kpc}$ and $\Lambda=10^{-22}\Lambda_{-22}\,{\rm erg}\,{\rm cm}^{3}\,{\rm s}^{-1}$ (see Fig. 1). With the exception of $\mathcal{B}$, the properties of the self-similar and integrated solutions differ by order-unity factors. Specifically, the integrated solution satisfies $t_{\rm flow}\approx t_{\rm cool}$ at all radii (lower-left panel), comparable to the constant $t_{\rm cool}/t_{\rm flow}=1.15$ in the self-similar solution (eqn. 22). The order unity differences between the integrated and self-similar solution are mainly due to the ‘decrease-flat-decrease’ shape of $T_{\rm c}(r)$, which is approximated as a straight power-law in the self-similar solution. The roughly constant offset in $\mathcal{B}$ between the two solutions is a result of a constant offset in $\Phi$, since the two solutions differ in the shape of $\Phi(r)$ at large scales.

The radii which limit the range of applicability of the solution, as discussed in section 2.2, are marked in the top-right panel. These radii are $R_{1/2}=3\,{\rm kpc}$ (eqn. 11), and $R_{\rm cool}=200\,{\rm kpc}$ derived from the condition $t_{\rm cool}=t_{\rm H}=13.6\,{\rm Gyr}$ ($t_{\rm H}$ is marked in the $t_{\rm cool}$ panel). We expect the derived solution not to be valid near or outside these limiting radii. Near $R_{1/2}$ a solution which includes angular-momentum must be found. Around $R_{\rm cool}$ one must find a time-dependent solution which accounts for the growth of $R_{\rm cool}$ with time. Bertschinger (1989) showed that as long as ${\rm d}R_{\rm cool}/{\rm d}t\ll c_{\rm s}(R_{\rm cool})$, such time-dependent solutions join smoothly onto the self-similar solutions. This property is demonstrated in section 3 using hydrodynamic simulations. If $R_{\rm shock}<R_{\rm cool}$ the maximum radius of applicability will be $R_{\rm shock}$, and at $R_{\rm shock}$ the solution must satisfy the shock jump conditions. In appendix B we demonstrate that such a solution is essentially identical to the marginally-bound solution within $R_{\rm shock}$.

2.5 Cooling flows in galaxy-scale halos

In this section we discuss several properties of cooling flows in halos characteristic of \sim$$L^{*} galaxies ($M_{\rm halo}\sim 10^{12}\,{\rm M_{\odot}}$), using both analytic estimates based on the self-similar solutions (section 2.3) and the more accurate integrated solutions (section 2.4). Figure 3 plots four integrated and four self-similar solutions for $Z=0.3\,{\rm Z_{\odot}}$, and $M_{12}=1$, where we define $M_{\rm halo}=10^{12}M_{12}\,{\rm M_{\odot}}$. For the integrated solutions we assume $R_{\rm sonic}=0.1\,{\rm kpc}$ (purple), $1\,{\rm kpc}$ (blue), $10\,{\rm kpc}$ (green), and $100\,{\rm kpc}$ (light green), while for the self-similar solutions (thin grey lines) we use the values of ${\dot{M}}$ found in the integrated solutions, as noted in the bottom-left panel. The blue solution with $R_{\rm sonic}=1\,{\rm kpc}$ is the solution shown in Fig. 2. The curves span the maximum range of applicability of the solutions, between $R_{1/2}$ and $R_{\rm cool}$. In cases where $R_{\rm sonic}>R_{1/2}$, dotted lines plot the inner supersonic part of the solution. We first address the subsonic part of the solutions which are the focus of this paper, and then address the supersonic part.

The self-similar solutions provide a good approximation to the subsonic part of the integrated solutions. Specifically, the integrated solutions satisfy $T\approx T_{\rm c}$ (top-left panel) and $t_{\rm flow}\approx t_{\rm cool}$ (bottom-middle), as expected from eqns. (21)–(22) derived in the context of the self-similar solutions. These two conditions demonstrate why cooling flows form a single-parameter family of solutions. The free parameter sets the normalization of other plotted properties – $n_{\rm H}$, $\mathcal{M}$, $t_{\rm cool}$, and $t_{\rm cool}/t_{\rm ff}$. The slope of the profiles of these properties is however roughly the same in all solutions.

The top-middle panel of Fig. 3 shows the density profile of the solutions. For the self-similar solutions, this profile can also be derived analytically. Approximating $v_{\rm c}$ as $140\,M_{12}^{0.36}r_{100}^{m}\,\rm km\ s^{-1}$ with $m=-0.1$ (Fig. 1), eqns. (21)–(22) imply $A=1.08$ and $B=0.87$, so we get from eqn. (27) that

[TABLE]

where ${\dot{M}}=1{\dot{M}}_{1}\,{\rm M_{\odot}}\,{\rm yr}^{-1}$, and we used $X=0.75$. This normalization of ${\dot{M}}$ is based on the SFR of the Milky Way (Bland-Hawthorn & Gerhard 2016), and corresponds to the purple lines in Fig. 3. Eqn. (29) implies that ${\dot{M}}\propto n_{\rm H}^{2}$, as expected in cooling flows since ${\dot{M}}\propto\rho v$ and $v\approx r/t_{\rm cool}\propto\rho$. Also, the scaling $n_{\rm H}\propto r^{-1.6}$ is consistent with estimates of the hot gas density slope in the Milky Way halo based on X-ray emission and absorption (Bregman et al. 2018). We return to this comparison in section 4.

The dashed black line in the top-middle panel plots the expected density if baryons follow the dark matter with a cosmic baryon fraction $f_{\rm b}$. This estimate is higher than the ${\dot{M}}=1\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ solution by an order of magnitude at $R_{\rm cool}$ (the right end of the plotted line), and by a larger factor at smaller radii. Hence the ${\dot{M}}=1\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ solution corresponds to a highly baryon-depleted halo.

Using eqns. (5), (24) and (29), we can derive an expression for $t_{\rm cool}$:

[TABLE]

which can be compared to the integrated solutions in the lower-left panel of Fig. 3. The cooling radius where $t_{\rm cool}=13.6\,{\rm Gyr}$ is hence

[TABLE]

Equation (31) provides an estimate of the outer limit of applicability of the cooling flow solution.

Similarly, the Mach number profile in cooling flows can be derived from eqn. (26):

[TABLE]

and the ratio $t_{\rm cool}/t_{\rm ff}$ in cooling flows follows from eqn. (28)

[TABLE]

Eqns. (32)–(33) imply that $t_{\rm cool}/t_{\rm ff}\propto\mathcal{M}^{-1}$ is a weak function of radius in \sim$$L^{*} halos, as can be seen in the right panels of Fig. 3. The flatness of $t_{\rm cool}/t_{\rm ff}$ is similar to the basic ansatz of thermal instability and ‘precipitation’ models based on simulations with galaxy feedback (Sharma et al. 2012b, Voit et al. 2017), though note that the physics are different since the cooling flow solution does not include feedback. We compare cooling flow solutions to precipitation models in section 5.

We now turn to discuss the supersonic parts of the solutions. A rough estimate of $R_{\rm sonic}$ where $\mathcal{M}=1$ can be derived from eqn. (32), which yields:

[TABLE]

This estimate can be compared to the integrated results in the top-right panel of Fig. 3. Within $R_{\rm sonic}$, the solutions rapidly lose their thermal energy (top-left panel) and the solution becomes a free-falling solution rather than a cooling flow solution. Within $R_{\rm sonic}$ the solution also has $t_{\rm cool}<t_{\rm ff}$ (bottom-right panel and eqn. 28), which is incompatible with the assumption of steady-state used to derive these solutions. That is, within $R_{\rm sonic}$ the flow is likely to be subject to thermal instabilities, repressurizing shocks, and inefficient mixing, so we do not expect the steady-state solutions to apply.

For solutions with $R_{\rm sonic}>R_{1/2}$ (e.g. the green and light green solutions in Fig. 3), one may ask whether the subsonic part of the solution at $R>R_{\rm sonic}$ is valid, if there is no valid steady-state solution within $R_{\rm sonic}$ which can support its weight against gravity? In Paper II we show that indeed, if $R_{\rm sonic}>R_{1/2}$ all the halo gas collapses on a dynamical timescale, including the subsonic gas beyond $R_{\rm sonic}$. Steady-state solutions for the halo gas are thus possible only if the gas remains subsonic down to $R_{1/2}$.

2.6 Cooling flows in cluster-scale halos

Figure 4 plots three cooling flow solutions in $10^{15}\,{\rm M_{\odot}}$ halos, derived assuming $Z=0.3\,{\rm Z_{\odot}}$ and either $R_{\rm sonic}=0.1\,{\rm kpc}$ (purple), $0.3\,{\rm kpc}$ (turquoise) or $1\,{\rm kpc}$ (green). The implied ${\dot{M}}$ are noted in the bottom-left panel. In cluster-scale halos the density at $R_{\rm cool}$ is directly observed via its X-ray emission (see section 4), and is found to be roughly consistent with the cosmic baryon fraction, corresponding to the normalization of the green solution in Fig. 4.

To derive the self-similar solution for cluster-scale halos, we approximate $v_{\rm c}$ as $850\,M_{15}^{0.23}r_{100}^{m}\,\rm km\ s^{-1}$ with $m\approx 0.3$, where $M_{\rm halo}=10^{15}M_{15}\,{\rm M_{\odot}}$ (see Fig. 1). In this case eqns. (21)–(22) imply $A=0.36$ and $B=1.4$. Hence in the self-similar solution we get (eqn. 24)

[TABLE]

The top-left panel in Fig. 4 shows that in the integrated solution $T$ is a factor of $1.5-2$ above $T_{\rm c}$, less than suggested by eqn. (35), and has a weaker dependence on radius than $T_{\rm c}$. These differences between the integrated and self-similar solutions occur because of the flattening of the $v_{\rm c}$ profile beyond $300\,{\rm kpc}$ (see Fig. 1), which affects the integrated solution but is not captured by the power-law approximation of the gravitational potential in the self-similar solution.

Calculating eqn. (27) in the context of clusters gives

[TABLE]

where we defined $\Lambda=10^{-23}\Lambda_{-23}\,\rm erg\,{\rm cm}^{3}\,{\rm s}^{-1}$ and ${\dot{M}}=10^{3}{\dot{M}}_{1000}\,{\rm M_{\odot}}\,{\rm yr}^{-1}$. For comparison, the integrated solutions suggest a somewhat steeper slope of $n_{\rm H}\propto r^{-1.4}$ (top-right panel of Fig. 4). This difference is also due to the flattening of $v_{\rm c}$ beyond $300\,{\rm kpc}$. Note that the density slope in cooling flows is similar to the dark matter slope at these radii. The cooling time is

[TABLE]

while the cooling radius is

[TABLE]

The value of $R_{\rm cool}$ is similar to that in \sim$$L^{*} galaxies (eqn. 30).

The entropy of the hot gas can be estimated from eqns. (35) and (36):

[TABLE]

where $n_{\rm e}\approx 1.2\,n_{\rm H}$ is the electron density. The derived slope is consistent with the slope found in the integrated solutions, shown in the bottom-left panel of Fig. 4. Moreover, this slope is steeper than the $K\propto r^{1.1}$ expected from self-similar cosmological accretion theory in the absence of cooling (e.g. Tozzi & Norman 2001). Therefore, in the absence of physical processes other than gravity and cooling, the entropy profile is expected to steepen around $R_{\rm cool}$ from $\sim r^{1.1}$ to $\sim r^{1.4}$. Note that the entropy profile scales as $\sim r^{1+4/3m}$ (eqn. 25) and $m$ becomes smaller with decreasing mass (Fig. 1), so the expected steepening at $R_{\rm cool}$ is weaker in group-scale halos than in cluster-scale halos.

Last, we calculate $t_{\rm cool}/t_{\rm ff}$ in clusters. Using eqn. (28) we get

[TABLE]

This relation can be compared to the bottom-right panel of Fig. 4, where we plot $t_{\rm cool}/t_{\rm ff}$ in the integrated solutions. Eqn. (40) implies that cooling flows in clusters are expected to remain stable ($t_{\rm cool}>t_{\rm ff}$) down to the galaxy scale $R_{1/2}\approx 30\,{\rm kpc}$.

3 Hydrodynamic simulations

In this section we present a series of idealized 3D hydrodynamic simulations of gas in dark matter halos. Concisely put, our simulations are controlled numerical experiments in which gas that is initially in hydrostatic equilibrium with an external gravitational potential is allowed to cool radiatively. Roughly hydrostatic conditions are expected in realistic halos after a virial shock has formed (see further discussion in Paper II). Our goal is to verify that in a hydrodynamic setting the halo gas indeed chooses one of the single-parameter family of solutions discussed in the previous section. The simulations are based on the simulations described in Fielding et al. (2017), though they do not include feedback from the central galaxy.

3.1 Simulation setup

The simulations are performed using the grid-based hydrodynamics code athena++ (Stone et al. in prep). We adopt an adiabatic equation of state with $\gamma=5/3$ and solve the standard hydrodynamics equations with additional source terms to include gravity and optically thin cooling and photoionization heating. We impose a cooling time constraint on the time step so that the global time step is the smaller of $\delta t_{\rm hydro}$ and ${\rm min}\{t_{\rm cool}\}/10$, to ensure that cells do not over-cool in one time step. Both the cooling and hydrodynamics updates are done with this time step. The same Wiersma et al. (2009) cooling and heating tables are used in the simulations as in the analytic calculations in the previous section (examples are shown in right panel of Fig. 1). Likewise, the simulations use the same gravitational accelerations that give rise to the $v_{\rm c}$ profiles shown in the left panel of Fig. 1.

3.1.1 Simulation domain and geometry

We adopt a spherical-polar coordinate system for our simulations. Given that our initial and boundary conditions are spherically-symmetric, we do not expect deviations from spherical symmetry to span broad angular ranges, so we restrict our computational domain to $\pi/4\leq\theta\leq 3\pi/4$ and $\pi/4\leq\phi\leq 3\pi/4$, where $\theta$ is the polar angle and $\phi$ is the azimuthal angle. This enables us to reduce the computation expense of a given simulation and to avoid the small time steps that arise from the Courant-Friedrichs-Lewy (CFL) condition near the poles. In the radial direction our domain extends from $3-1000\,{\rm kpc}$ with logarithmic grid spacing. The outer radial boundary is chosen to be beyond $R_{\rm cool}$ (eqns. 31 and 38), while the inner boundary is chosen to be within the galaxy scale for $\gtrsim 10^{12}\,{\rm M_{\odot}}$ halos (eqn. 11). Our fiducial simulations have approximately 1:1 cell aspect ratios with $64$ cells in the angular directions and $240$ in the radial direction. This corresponds to cell widths ranging from $73\,{\rm pc}$ at the inner boundary, $720\,{\rm pc}$ at $R=30\,{\rm kpc}$, and $24\,{\rm kpc}$ at the outer boundary. Simulations with twice and half the resolution resulted in nearly indistinguishable radial profiles of the gas properties, indicating that the simulations are well converged at our fiducial resolution.

3.1.2 Boundary conditions

A standard radial outflow boundary condition, which imposes zero gradients, is insufficient for our problem because the cooling flow solution that develops has non-zero gradients. Thus, if the flow is subsonic at the inner boundary (i.e. if $R_{\rm sonic}<R_{\rm inner}=3\,{\rm kpc}$) then a zero gradient boundary precludes the development of a steady-state solution. We therefore adopt boundary conditions that quadratically extrapolates the primitive variables into the ghost zones. In the polar and azimuthal directions we adopt periodic boundary conditions.

In practice, the inner radial boundary condition acts as a sink that allows gas to flow through it without introducing significant numerical artifacts. The outer boundary is beyond $R_{\rm cool}$ in all cases, so there are few changes near the outer boundary over the durations of the simulations, and the extrapolation boundary condition primarily serves to maintain equilibrium with the gravitational potential.

3.1.3 Initial conditions

We run simulations with a $z=0$ UV background and a range of halo masses, $M_{\rm halo}$, and use the appropriate $v_{\rm c}(r)$ plotted in Fig. 1. We start with an initially static halo ($v(r)=0$), with the density and temperature given by the hydrostatic relation

[TABLE]

One needs also to specify the entropy profile and the total gas mass. As any bound hydrostatic solution has $c_{\rm s}^{2}\sim v_{\rm c}^{2}$, for simplicity we assume $c_{\rm s}^{2}=v_{\rm c}^{2}$, which yields an entropy profile $K\propto r^{\gamma(\gamma-1)}\sim r^{1.1}$ at radii where $v_{\rm c}(r)$ is flat. We checked that our conclusions are not sensitive to the exact choice of the ratio $v_{\rm c}^{2}/c_{\rm s}^{2}$, as long as it is of order unity. The total gas mass within $R_{\rm vir}$ is set to equal some fraction of baryon budget $f_{\rm b}M_{\rm halo}$, as detailed below.

To verify that multiphase structure does not develop significantly in the simulation, as suggested by linear theory when $t_{\rm cool}>t_{\rm ff}$ (e.g. Balbus & Soker 1989), we impose isobaric density perturbations in the initial conditions. The perturbations are generated so that the power is evenly distributed ($P_{\delta\rho}\propto k^{0}$) between modes satisfying $1\leq kR_{\rm outer}/2\pi\leq 64$ with an amplitude so that $\langle\delta\rho/\rho\rangle_{\rm rms}=0.1$, where $R_{\rm outer}$ is the outer radius of the domain. Testing indicates that our results are not sensitive to the exact choices we made for the perturbations, except where noted below.

3.2 Results

Figure 5 plots the results of four different simulations, with $M_{\rm halo}=10^{15}\,{\rm M_{\odot}}$ (left column), $10^{13}\,{\rm M_{\odot}}$ (second column) and $10^{12}\,{\rm M_{\odot}}$ (two right columns). The first three simulations initially have $M_{\rm gas}=f_{\rm b}M_{\rm halo}-M_{*}$ within $R_{\rm vir}$, while the simulation on the right starts with a gas mass equal to a third of this value. A metallicity of $\,{\rm Z_{\odot}}/3$ is assumed throughout. In each column we show the temperature in units of $T_{\rm vir}$ (eqn. 9) and the radial velocity, at different simulation times $t$ as noted in the panels. To project the snapshot onto the plotted $r-\theta$ plane we calculate the mass-weighted average of $T$ and $v$ over the simulated range in $\phi$. The initial hydrostatic conditions are plotted in the top row. Note how an inflow develops in all simulations. Even though we seed perturbations in all simulations, strong deviations from spherical symmetry are apparent only at late times in the fiducial density $10^{12}\,{\rm M_{\odot}}$ simulation, and near the inner boundary. As we show below, only at these regions and times the gas is supersonic and has $t_{\rm cool}<t_{\rm ff}$, so multi-phase structure can develop.

Solid colored lines in Figure 6 plot shell-averaged properties in the $10^{15}\,{\rm M_{\odot}}$ simulation, as a function of radius and time. The panels show $T$ (top-left), $n_{\rm H}$ (top-right), $\mathcal{M}$ (middle-left), ${\dot{M}}$ (middle-right), $t_{\rm cool}/t_{\rm ff}$ (bottom-left), and $t_{\rm cool}/t_{\rm flow}$ (bottom-right). For each shell with radius $R$, we calculate mass-weighted averages for $T$, $n_{\rm H}$, $\mathcal{M}$, and $v$, while ${\dot{M}}$ is the volume integral of $\rho v$ in the shell. The shell cooling time is calculated as

[TABLE]

where ${\rm d}\Omega=\sin\theta{\rm d}\theta{\rm d}\phi$. Eqn. (42) implies that the shell cooling time is the ratio of the total thermal energy in the shell to the total luminosity of the shell.

Vertical ticks in Fig. 6 mark the cooling radius $R_{\rm cool}$ at the different snapshots, defined via $t_{\rm cool}(R_{\rm cool},t)=t$. This radius mildly increases from $R_{\rm cool}=50\,{\rm kpc}$ at $t=1\,{\rm Gyr}$ to $R_{\rm cool}=130\,{\rm kpc}$ at $t=10\,{\rm Gyr}$. The top two panels demonstrate that within $R_{\rm cool}$, $T$ increases and $n_{\rm H}$ decreases relative to their initial values (light green lines), while they are essentially unchanged at larger radii.

The two dashed lines plot steady-state solutions with $R_{\rm sonic}=0.6\,{\rm kpc}$ (blue) and $0.3\,{\rm kpc}$ (purple), derived as discussed in section 2. The values of $R_{\rm sonic}$ are chosen so ${\dot{M}}$ in the steady-state solution equals ${\dot{M}}$ in the $t=3\,{\rm Gyr}$ and $t=10\,{\rm Gyr}$ snapshots, respectively. Note that only the subsonic part of these two transonic solutions is within the simulated domain. The radial profiles of all plotted properties are well reproduced by the steady-state solutions, between $\approx 6\,{\rm kpc}$ and $\approx R_{\rm cool}$. This result supports the main point of this paper, that cooling gaseous halos converge onto the steady-state cooling flow solutions described above, and that the entire solution is determined with a single free parameter.

The deviation of the radial profiles from the steady-state solutions within $\approx 6\,{\rm kpc}$ is a boundary effect, since our boundary conditions are not exactly consistent with a steady-state solution (see section 3.1.2). Other choices of the inner boundary radius support this conclusion – the radial profiles in the simulation deviate from a steady-state solution within a factor of $\approx 2$ from the inner boundary.

The lower-right panel in Fig. 6 shows that $t_{\rm cool}\approx t_{\rm flow}$ out to the outer boundary of the simulation, including radii larger than $R_{\rm cool}$ where $t_{\rm cool}\gg t$. This behavior is somewhat surprising given that a cooling flow is expected to form on a timescale $t_{\rm cool}$. In Appendix C we demonstrate why the relation $t_{\rm cool}\approx t_{\rm flow}$ is established on a sound crossing timescale rather than a cooling timescale. The condition $t_{\rm cool}\approx t_{\rm flow}$ is though insufficient for a flow to lie on one of the steady-state solutions. The steady-state solutions also require a constant ${\dot{M}}$ and the appropriate entropy profile (eqns. 1 and 3), which are both established on a timescale $t_{\rm cool}$.

What sets the free parameter ${\dot{M}}$ in the simulation, which in Fig. 6 decreases from $520\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ at $t=3\,{\rm Gyr}$ to $240\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ at $t=10\,{\rm Gyr}$? Since the conditions beyond $R_{\rm cool}$ remain near their initial values, the free parameter must allow the cooling flow solution within $R_{\rm cool}$ to smoothly join the roughly static initial conditions beyond $R_{\rm cool}$. The value of ${\dot{M}}(t)$ is hence set by the initial conditions at a radius $R_{\rm cool}(t)$. This relation was worked out analytically by Bertschinger (1989) for the self-similar case, who found that up to a factor of order unity

[TABLE]

where $\rho_{0}(R)$ is the initial gas density at radius $R$. We find that ${\dot{M}}$ near the inner boundary (at $R=20\,{\rm kpc}$) in the $10^{15}\,{\rm M_{\odot}}$ simulation is consistent with eqn. (43) to within a factor of $1.5$. For comparison, the more naive estimate of ${\dot{M}}(t)=M_{\rm gas}(<R_{\rm cool}(t))/t$ is larger by a factor of $\sim 3$ then suggested by eqn. (43). Note that if the flow is bounded on the outside by an accretion shock rather than by $R_{\rm cool}$ (see section 2.2), the jump conditions at the shock will determine ${\dot{M}}(t)$.

Figures 7 and 8 plot the spherically-averaged properties of the $10^{13}\,{\rm M_{\odot}}$ simulation and the low-density $10^{12}\,{\rm M_{\odot}}$ simulation, together with steady-state solutions corresponding to the $t=3\,{\rm Gyr}$ and $t=10\,{\rm Gyr}$ snapshots. As in Fig. 6, the spherically-averaged properties of the simulation are well reproduced by the steady-state solutions, between $\approx 6\,{\rm kpc}$ and out to $\approx R_{\rm cool}$. These simulations hence also support the main conclusion of this study, that within $R_{\rm cool}$ the halo gas settles on one of the single-parameter solutions derived in section 2. The value of ${\dot{M}}(t)$ near the inner boundary in these simulations is consistent with eqn. (43) to within a factor of two.

Figure 9 plots average radial profiles in the $10^{12}\,{\rm M_{\odot}}$ simulation with an initially closed baryon fraction (third column in Fig. 5). To increase the dynamic range we use a ‘symmetric log’ axis in the $\mathcal{M}$ and ${\dot{M}}$ panels, where the y-axis is linear for absolute values smaller than unity and logarithmic for absolute values above unity. At $t=3\,{\rm Gyr}$, the flow properties are consistent with the steady-state solution in the range $6\,{\rm kpc}<R<R_{\rm cool}$, as seen in the snapshots plotted in Figs. 6 – 8. However, at $t=10\,{\rm Gyr}$ the flow is supersonic, with $t_{\rm cool}<t_{\rm ff}$ and a temperature significantly below virial. Both $R_{\rm sonic}$ and $R_{\rm cool}$ in this snapshot are at $800\,{\rm kpc}$. The steady-state solutions derived in the previous section are thus invalid if $t_{\rm cool}<t_{\rm ff}$. The transition between a cooling flow and a supersonic inflow as seen in this simulation is discussed in Paper II.

Figure 10 plots the r.m.s. density dispersion $\langle\delta\rho/\rho\rangle_{\rm rms}$ in radial shells in the four simulations. Each panel corresponds to a different simulation, while each colored line corresponds to a different snapshot as noted in the legend in the top-right panel. The profiles of $t_{\rm ff}/t_{\rm cool}$ in the $t=10\,{\rm Gyr}$ snapshots are plotted with dashed grey lines. In the $10^{15}\,{\rm M_{\odot}}$ simulation (top-left panel), by $t=1\,{\rm Gyr}$ the amplitude of the perturbations has decreased from their inital value of $\langle\delta\rho/\rho\rangle_{\rm rms}=0.1$ (light green lines). Specifically, at $t=1$, $3$ and $10\,{\rm Gyr}$ the relative fluctuations are significantly below unity, indicating that multi-phase structure has not developed. A similar behavior is apparent in the $10^{13}\,{\rm M_{\odot}}$ and low density $10^{12}\,{\rm M_{\odot}}$ simulations (right panels) except at small $t$ and large $R$ where the initial perturbations remain. Significant fluctuations where $\langle\delta\rho/\rho\rangle_{\rm rms}$ is larger than unity are apparent only in the fiducial density $10^{12}\,{\rm M_{\odot}}$ simulation at $t=10\,{\rm Gyr}$, in which $t_{\rm ff}/t_{\rm cool}>1$ and the flow is supersonic. These fluctuations are also apparent in the corresponding temperature panel in Fig. 5. Figure 10 thus demonstrates that cooling flows do not develop into a multiphase CGM, as long as $t_{\rm cool}>t_{\rm ff}$.

We emphasize that the latter conclusion applies only to perturbations with an initial amplitude significantly lower than unity. Stronger perturbations seeded e.g. by the wakes of satellite galaxies are likely to persist in the halo, and may affect the evolution of the ambient medium even if it has $t_{\rm cool}\gtrsim t_{\rm ff}$ (e.g. Sharma et al. 2012a; Choudhury et al. 2019).

The decrease in amplitude of the perturbations from their initial value occurs roughly on a free-fall time, which increases with decreasing halo mass at fixed $r$. After this initial phase, in the three simulations which do not collapse the amplitude of the perturbations is within a factor of $2-3$ of $t_{\rm ff}/t_{\rm cool}$ (compare the purple lines with the dashed gray lines). We find an almost identical final perturbation amplitude in simulations initialized with a higher $\langle\delta\rho/\rho\rangle_{\rm rms}=0.3$, though simulations with a lower initial $\langle\delta\rho/\rho\rangle_{\rm rms}=0.01$ have a final $\langle\delta\rho/\rho\rangle_{\rm rms}\ll t_{\rm ff}/t_{\rm cool}$. A similar saturation of density fluctuations at an amplitude $\approx t_{\rm ff}/t_{\rm cool}$ has previously been seen in simulations of thermally balanced atmospheres, and attributed to dissipation of perturbations of this amplitude via non-linear mode coupling (McCourt et al. 2012; Voit et al. 2017).

4 Comparison with Observations

In this section we compare observables of cooling flow solutions with observational constraints of halo gas at low redshift.

4.1 OVII and OVIII absorption in the Milky Way

The expected column of an ion $X^{i+}$ in the halo can be derived via

[TABLE]

where $N_{X^{i+}}$ is the ion column, ${\rm X}/{\rm H}$ is the abundance of element ${\rm X}$ relative to hydrogen, $f_{\rm X^{i+}}$ is the fraction of X particles in the $i$-th ionization state, and ${\rm d}s$ is the line-of-sight element. For $f_{\rm X^{i+}}$ we use the ionization equilibrium calculations in the trident package (Hummels et al. 2017), which used cloudy (Ferland et al. 2013) and assumed the gas is exposed to a HM12 background. As in section 2.1, we assume the actual background has twice the intensity deduced by HM12.

We calculate two cooling flow solutions with a metallicity of either $0.3\,{\rm Z_{\odot}}$ or $\,{\rm Z_{\odot}}$, ${\dot{M}}=1.65\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ equal to the SFR of the Galaxy, $M_{\rm halo}=1.3\cdot 10^{12}\,{\rm M_{\odot}}$, and $M_{*}=6.3\cdot 10^{10}\,{\rm M_{\odot}}$. All estimates for Galactic parameters are taken from Bland-Hawthorn & Gerhard (2016). The choice of ${\dot{M}}={\rm SFR}$ is motivated by the assumption that all fuel for star formation is provided by the hot CGM, though lower ${\dot{M}}$ are also possible in the context of a cooling flow CGM, if some of the fuel is supplied by stellar mass loss (Leitner & Kravtsov 2011). The values of the metallicity are chosen to bracket the possible range, where the lower limit was deduced by Miller & Bregman (2013, 2015), who divided the observed $N_{\text{\sc O\,vii}}$ by the upper limit on the dispersion measure towards the LMC (see section 4.3 below), while the upper limit is the ISM metallicity.

Using $T(r)$ and $n_{\rm H}(r)$ from the cooling flow solutions (where $r$ is the galactocentric radius), we integrate eqn. (44) along different sightlines through the MW halo. For a given solution, the predicted columns depend on the angle of the sightline relative to the Galactic center. We end the integration at $R_{\rm cool}=130\,{\rm kpc}$ (for the $0.3\,{\rm Z_{\odot}}$ solution) or at $R_{\rm cool}=180\,{\rm kpc}$ ($\,{\rm Z_{\odot}}$), though since more than 85% of the absorption originates within the inner $\sim 40\,{\rm kpc}$, the outer end of the integration does not significantly affect the predicted columns. The predicted $N_{\text{\sc O\,vii}}$ and $N_{\text{\sc O\,viii}}$ are shown in the top-left and top-right panels of Figure 11, respectively.

The observed $N_{\rm OVII}$ plotted in the top-left panel of Fig. 11 are taken from Gupta et al. (2012, hereafter G12), who presented a sample of eight high S/N sightlines with detections in both the K$\alpha$ (21.6Å) and K$\beta$ (18.63Å) absorption lines. For comparison, Faerman et al. (2017) deduced a similar median $N_{\text{\sc O\,vii}}=1.4\cdot 10^{16}\,{\rm cm}^{-2}$ based on the sample of Fang et al. (2015), which included 43 sightlines with lower S/N than G12. The $N_{\rm OVII}$ measurements in Fang et al. (2015) are based solely on the K$\alpha$ line and typically have \gtrsim$$1\,{\rm dex} uncertainties, so we do not show them individually. The bottom-left panel shows the difference between the observations and the predictions of the cooling flow solutions. The cooling flow solution with $Z=\,{\rm Z_{\odot}}/3$ is consistent with the G12 observations (reduced $\chi^{2}=0.8$) while the $Z=\,{\rm Z_{\odot}}$ solution is marginally consistent (reduced $\chi^{2}=1.7$). All solutions with metallicities $0.3-0.7\,{\rm Z_{\odot}}$ yield an acceptable fit with a reduced $\chi^{2}$ lower than unity.

The observed $N_{\rm OVIII}$ shown in the top-right panel of Fig. 11 are also based on the G12 data. We calculate $N_{\rm OVIII}$ using a curve-of-growth analysis, based on the O viii equivalent widths and the velocity width parameter $b$ measured for O vii, which is justified if the two ions originate in the same gas (see similar calculation in Faerman et al. 2017). Comparison of these measurements with the models yields a reduced $\chi^{2}=0.7$ for the $\,{\rm Z_{\odot}}$ model and a reduced $\chi^{2}=1.1$ for the $\,{\rm Z_{\odot}}/3$ model. Combining the $N_{\rm OVII}$ and $N_{\rm OVIII}$ observations, the cooling flow solutions yield acceptable fits for metallicities in the range $0.3-0.8\,{\rm Z_{\odot}}$. Fig. 11 thus demonstrates that the cooling flow solutions are consistent with the observed $N_{\rm OVII}$ and $N_{\rm OVIII}$. We emphasize that there are no free parameters in these solutions beyond the uncertainty in $Z$.

The cooling flow solutions used to calculate the predicted columns in Fig. 11 assume a mass inflow rate equal to the SFR in the Milky Way. We also calculate solutions with different values of ${\dot{M}}$ and compare them to the observational data. We find an acceptable fit for $1.6<{\dot{M}}<4\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ if $Z=\,{\rm Z_{\odot}}/3$, and for $0.2<{\dot{M}}<1.3\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ if $Z=\,{\rm Z_{\odot}}$. Therefore, if the MW halo gas forms a cooling flow, then the observed X-ray ion columns allow for mass inflow rates in the range $0.2-4\,{\rm M_{\odot}}\,{\rm yr}^{-1}$, or equivalently $0.12-2.4$ times the SFR.

4.2 OVII and OVIII emission in the Milky Way

Li & Bregman (2017, hereafter LB17) derived the density profile of the Milky Way halo gas based on X-ray observations of O vii and O viii emission lines, using 648 sightlines from the sample of Henley & Shelton (2012, 2013). In order to derive the density profile they assumed a constant temperature of $2\cdot 10^{6}\,{\rm K}$, similar to the cooling flow solutions for the MW halo discussed in the previous section, in which the temperature can be approximated as $T/10^{6}\,{\rm K}\approx 1.6-0.5\log(R/10\,{\rm kpc})$. The observationally-constrained density profile of LB17 can thus be directly compared to the density profile predicted by cooling flow solutions.

The most detailed model in LB17 (model ‘9’), which accounted for radiation transfer effects, included a disk component and allowed for rotation of the hot gas, yielded an electron density profile at large radii of

[TABLE]

This density profile was constrained using the O vii emission line triplet near 22Å. LB17 also constrained their model based on the O vii and O viii emission lines near 19Å (their model ‘3’), which yielded a profile consistent with eqn. (45), though with larger errors. The top panel of Figure 12 plots the LB17 density profile versus the MW cooling flow solution with ${\dot{M}}={\rm SFR}$ and $Z=0.3\,{\rm Z_{\odot}}$. LB17 also assumed $Z=0.3\,{\rm Z_{\odot}}$, though since the density in cooling flows and in LB17 scale similarly with the assumed $Z$, the assumed $Z$ does not affect the comparison. Fig. 12 shows that the slope of the LB17 model and the cooling flow solutions are almost identical at radii $\gtrsim 10\,{\rm kpc}$, as also evident from comparing eqn. (45) with eqn. (29). As the density slope is a robust prediction of the cooling flows, the LB17 result supports the hypothesis that the MW halo gas beyond $\gtrsim 10\,{\rm kpc}$ forms a cooling flow. At smaller radii the LB17 profile flattens while the cooling flow profile remains steep, though we do not expect the cooling flow to apply at such small radii due to angular momentum and deviations of the potential from spherical symmetry.

Fig. 12 shows that the normalization of the LB17 and cooling flow profiles differ by a factor of $\approx 2$. Roughly $30\%$ of this difference can be accounted for by the higher temperature assumed in LB17 compared to the temperatures in the cooling flow solution, which yields a higher emissivity and hence a lower density for a given line emission. The remaining difference may suggest that the actual mass inflow rate is a factor of $\approx 2$ lower than the value of ${\rm SFR}=1.65\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ assumed in the cooling flow solution.

To produce a good fit to the observations, LB17 assumed an uncertainty of 2.1 L.U. in their model for the O vii emission ($\approx 40\%$ of the median O vii intensity). They attributed this uncertainty to unaccounted variations in the emission from either the Local Bubble or the halo gas. If these emission variations are not due to the Local Bubble but arise in halo gas, they could indicate density variations of order $\lesssim 20\%$ in the cooling flow relative to spherical symmetry (since emission scales as $n_{\rm H}^{2}$).

In the bottom panel of Fig. 12 we show the total hot gas mass in the MW halo implied by the cooling flow solutions. The solutions are extended beyond $R_{\rm cool}$ (marked by a tick) assuming a hydrostatic pressure profile with $K\propto R$, though we note that this extension is uncertain. The implied total gas mass within $R_{\rm vir}$ is $17\%$ of the halo baryon budget. Extending the profile beyond $R_{\rm cool}$ with a shallower entropy profile yields a somewhat higher gas mass (e.g. $K\propto R^{0.5}$ yields a gas mass of $22\%$ of the baryon budget), while assuming a metallicity higher than $0.3\,{\rm Z_{\odot}}$ decreases the implied gas mass by $\approx(Z/0.3\,{\rm Z_{\odot}})^{-1/2}$. A similar low baryon mass was derived by Bregman et al. (2018), who deduced the total hot gas mass by extrapolating the LB17 gas density profile out to $R_{\rm vir}$ (marked with diamonds in Fig. 12).

4.3 Dispersion measure towards the LMC

Using measurements of the dispersion measure towards pulsars in the Large Magellanic Cloud (LMC), Anderson & Bregman (2010) estimated an electron column of $23\,{\rm cm}^{-3}\,{\rm pc}$ in the halo, after subtracting an estimated contribution of $47\,{\rm cm}^{-3}\,{\rm pc}$ from the Galactic disk and assuming a negligible contribution from gas in the LMC. For comparison, the cooling flow solution with ${\dot{M}}={\rm SFR}$ and $Z=0.3\,{\rm Z_{\odot}}$ predicts a dispersion measure between $8.3$ and $50\,{\rm kpc}$ of $37\,{\rm cm}^{-3}\,{\rm pc}$, while the $Z=\,{\rm Z_{\odot}}$ solution predicts $22\,{\rm cm}^{-3}\,{\rm pc}$. The predicted and observed values are comparable, supporting the cooling flow solution for hot gas in the MW halo. The somewhat lower values suggested by the dispersion measure observations, especially if $Z$ is relatively low or if the LMC contribution is non-negligible, may suggest ${\dot{M}}$ is actually somewhat smaller than the SFR, as also suggested by the comparison of the cooling flow solution with the LB17 profile in Fig. 12.

4.4 OVI absorption

Figure 13 compares the prediction of cooling flow solutions with observations of O vi absorption around $z\sim 0.2$ star-forming galaxies and around the Milky Way. Blue markers in the four left panels denote $N_{\text{\sc O\,vi}}$ measurements from Werk et al. (2013), Johnson et al. (2015), and Johnson et al. (2017). The galaxies are grouped by $M_{\rm halo}$ as noted at the top of each panel, where $M_{\rm halo}$ is estimated from the stellar masses of the central galaxies noted in the papers, using the Behroozi et al. (2018) relation111One object, J1435+3604_68_12, has $M_{\rm halo}=10^{13}\,{\rm M_{\odot}}$ and is not shown. O vi is not detected in this object.. For each $M_{\rm halo}$ bin, we calculate cooling flow models with the median $z=0.2$, either $Z=0.3\,{\rm Z_{\odot}}$ or $Z=\,{\rm Z_{\odot}}$, and ${\dot{M}}$ equal to the median SFR in the observed galaxies, taken from Werk et al. (2013). In the left panel where all galaxies do not have a published SFR estimate we assume ${\dot{M}}=0.07\,{\rm M_{\odot}}\,{\rm yr}^{-1}$, the average SFR for galaxies with the median $M_{*}=1.6\cdot 10^{8}\,{\rm M_{\odot}}$ in this bin (Behroozi et al. 2018). In the right panel the lower marker shows the average of observations along MW sightlines from Sembach et al. (2003), multiplied by two to mimic an external galaxy sightline as in the other panels. Due to confusion with absorption in the Galactic disk, the Sembach et al. (2003) observations do not include O vi absorption at local standard of rest velocities $\left|v_{\rm LSR}\right|<100\,\rm km\ s^{-1}$, which suggests $N_{\text{\sc O\,vi}}$ are actually a factor of $\sim 2$ higher (Zheng et al. 2015). The corrected average $N_{\text{\sc O\,vi}}$ is noted with the upper marker. As mentioned by Zheng et al. (2015), the corrected $N_{\text{\sc O\,vi}}$ in the Milky Way are similar to columns observed in other star-forming galaxies with the same mass, suggesting a similar physical origin. The cooling flow solutions in the right panel are the same solutions used in Fig. 11.

The predicted ion columns are calculated from the cooling flow solutions using eqn. (44), integrated along sightlines with different impact parameters $R_{\perp}$. We integrate out to $2R_{\rm vir}$ rather than to $R_{\rm cool}$ to avoid an unphysical break in the predicted $N_{\text{\sc O\,vi}}(R_{\perp})$ profile near $R_{\rm cool}$, though in all cases the contribution to $N_{\text{\sc O\,vi}}$ of gas beyond $R_{\rm cool}$ is small for sightlines with $R_{\perp}<R_{\rm cool}$. The cooling flows solutions underpredict the observed $N_{\text{\sc O\,vi}}$ around most galaxies, typically by a factor of $\sim 5$. A similar conclusion arises from the MW observations (right panel) where the solutions underpredict the average Sembach et al. (2003) observations by a factor of $1.8-3.2$, and the corrected observations from Zheng et al. (2015) by a factor of $3.5-6.5$. This failure in reproducing the observed $N_{\text{\sc O\,vi}}$ is in contrast with the cooling flow solutions success in reproducing the observed O vii and O viii absorption (Fig. 11), the density slope suggested by the observed O vii and O viii emission (Fig. 12), and the dispersion measure towards the LMC (section 4.3). In the discussion we suggest possible resolutions to this apparent discrepancy, in which either the O vi-bearing outer halo is preferentially heated, or alternatively the cooling flow extends only out to an accretion shock at $\lesssim 100\,{\rm kpc}$ and O vi predominantly traces cool gas beyond the shock.

The predicted $N_{\text{\sc O\,vi}}$ in cooling flows scales roughly as $N_{\text{\sc O\,vi}}\propto n_{\rm H}f_{\text{\sc O\,vi}}\propto{\dot{M}}^{1/2}$, since $n_{\rm H}\propto{\dot{M}}^{1/2}$ (eqn. 29) while $f_{\text{\sc O\,vi}}$ which depends mainly on $T$ is roughly independent of ${\dot{M}}$. Hence, in principle a cooling flow with ${\dot{M}}\sim 25\,{\rm SFR}$ could reproduce the O vi observations, though such a model would have $t_{\rm cool}<t_{\rm ff}$ and hence would be unstable. It has been shown that any model which assumes O vi traces radiatively cooling gas requires ${\dot{M}}\gg{\rm SFR}$ (Mathews & Prochaska 2017; Faerman et al. 2017; McQuinn & Werk 2018; Stern et al. 2018).

4.5 Cool clouds

Cool gas clouds ($\sim 10^{4}\,{\rm K}$) are routinely observed in dark matter halos via their low-ion (e.g. Mg ii, C ii, Si ii) and H i absorption (see review by Tumlinson et al. 2017). Around low-redshift \sim$$L^{*} galaxies, such absorption features appear in $\gtrsim 50\%$ of sightlines to background quasars with impact parameter out to $\approx 0.5R_{\rm vir}$ (Werk et al. 2013; Liang & Chen 2014, see also figure 11 in Stern et al. 2018). The prevalence of cool gas in halos is a challenge for the validity of subsonic cooling flow solutions, in which the halo gas is expected to be predominantly single-phased. This expectation is evident in Fig. 10, which demonstrates that in simulations seeded with small amplitude perturbations ($\langle\delta\rho/\rho\rangle_{\rm rms}=0.1$, see section 3.1.3) the density dispersion remains smaller than unity in all snapshots where $t_{\rm cool}>t_{\rm ff}$, consistent with results based on linear perturbation theory (Malagoli et al. 1987; Balbus & Soker 1989).

It is however important to note that due to projection effects, the large observed covering factor of cool clouds does not necessarily imply that they are widespread throughout the halo. Cool clouds clustered in a small fraction of the halo volume near ‘local’ disturbances which produce non-linear perturbations, such as in the extended disk, along collimated outflows, or near satellites and the material stripped from them (e.g. the Magellanic Stream), could in principle have a large area covering factor. This association of cool clouds with specific locations in the halo is suggested by the tendency of Mg ii clouds to align with either the minor axis or the major axis of the galaxy (Bouché et al. 2012; Kacprzak et al. 2012; Nielsen et al. 2015; Martin et al. 2019). It is thus a prediction of the cooling flow scenario that the observed cool clouds are limited to a small fraction of the halo volume. Such a scenario can be seen in the FIRE simulations where a large fraction of the cool gas in $\gtrsim 10^{12}\,{\rm M_{\odot}}$ halos is associated with satellite galaxies and their outflows (Faucher-Giguère et al. 2015, 2016), and in the simulations of Hummels et al. (2018), where cool clouds form only near the galaxy or near an inflowing filament.

In regions where cool clouds do form, they can be used as a rough barometer of the ambient virial-temperature gas by assuming pressure equilibrium between the cool and hot phases, as initially done by Spitzer (1956). In the top panel of Fig. 12 we plot the hot gas density implied by the pressure of clouds in the Magellanic stream derived by Stanimirović et al. (2002), who used 21cm observations and deduced $P/k=300\,{\rm cm}^{-3}\,{\rm K}$ assuming a distance of $45\,{\rm kpc}$. We convert their pressure estimate to a density estimate using $T(45\,{\rm kpc})=1.3\cdot 10^{6}\,{\rm K}$ in the cooling flow solution. We estimate an error of $0.3\,{\rm dex}$ on the Stanimirović et al. (2002) estimate based on the dispersion between different Magellanic stream clouds (see their figure 16). This pressure estimate is consistent with the cooling flow solution.

Also plotted in Fig. 12 is the hot gas density implied by the pressure in the local interstellar medium, as measured by Jenkins & Tripp (2011). They used observations of C i absorption and deduced a pressure of $P/k=3800\,{\rm cm}^{-3}\,{\rm K}$ with a dispersion of $0.175\,{\rm dex}$, which we convert to a density estimate using $T(8.3\,{\rm kpc})=1.6\cdot 10^{6}\,{\rm K}$ in the cooling flow solution. This pressure estimate is marginally consistent with the $0.3\,{\rm Z_{\odot}}$ cooling flow solution. The median value of Jenkins & Tripp (2011) is equal to the prediction of the $Z=\,{\rm Z_{\odot}}$ cooling flow solution (not shown).

The pressure of cool clouds in halos can also be constrained via photoionization modelling of low-ionization UV absorption features. In the COS-Halos sample of low-redshift star-forming galaxies, Werk et al. (2014) and McQuinn & Werk (2018) deduced a gas pressure which is a factor of $\gtrsim 10$ lower than expected in a gaseous halo with a closed baryon fraction. Broadly consistent results were obtained on a similar sample by Keeney et al. (2017, see also Voit et al. 2019). Though these estimates are subject to several systematic uncertainties (Stern et al. 2016; Chen et al. 2017), taking these pressure estimates at face value suggests a hot gas mass significantly below the baryon budget, roughly consistent with the $\lesssim 20\%$ baryon fraction deduced from the cooling flow solutions in section 4.2.

4.6 Clusters

The top panel in Figure 14 compares the predicted entropy profile in cooling flows with deprojected observations of five bright cool-core clusters from Hogan et al. (2017). The dashed cyan line plots the scaling $K\propto r^{1.1}$ expected in self-similar cosmological accretion, i.e. when cooling is neglected (Tozzi & Norman 2001; Voit et al. 2005). This prediction provides a good fit to the observations at large radii. The thick black line plots a cooling flow solution with $M_{\rm halo}=10^{15}\,{\rm M_{\odot}}$, $Z=0.3\,{\rm Z_{\odot}}$, and ${\dot{M}}=300\,{\rm M_{\odot}}\,{\rm yr}^{-1}$, where the latter is chosen to match the observations at $R_{\rm cool}=160\,{\rm kpc}$. The plotted curve spans between $R_{1/2}=30\,{\rm kpc}$ (eqn. 11) and $R_{\rm cool}$. The entropy profile in the cooling flow scales as $r^{1.4}$ (eqn. 39), steeper than expected from gravity alone. Hence, without additional physical processes beyond gravity and cooling, the entropy profile is expected to steepen at $R_{\rm cool}$. In contrast, the observed entropy profile appears to flatten, reaching $K\propto r^{0.7}$ at small scales (Panagoulia et al. 2014; Babyk et al. 2018). The flattening of the profile appears to occur near $R_{\rm cool}$, suggesting that feedback affects the intracluster medium out to the maximum radius where it is actively cooling.

The Hogan et al. (2017) clusters shown in the top panel were chosen mainly for their brightness, and are hence relatively representative of cool-core clusters. The bottom panel shows model fits to four objects from McDonald et al. (2019a), which were selected to be the most dynamically-relaxed clusters out of the 100 clusters in the South Pole Telescope sample (SPT, Bleem et al. 2015). The entropy profiles at $30-160\,{\rm kpc}$ of these clusters is steeper than in the sample shown in the top panel, and similar to the slope predicted by the cooling flow solution, especially in the Phoenix and SPT-CLJ2043-5035 clusters. The intracluster medium (ICM) in the most dynamically relaxed clusters may hence be forming a cooling flow at these radii, i.e. it may not be subject to significant heating by feedback. The entropy profiles flatten relative to the cooling flow solution within $\sim 30\,{\rm kpc}\sim R_{1/2}$.

In recent work published after the submission of this manuscript, McDonald et al. (2019b) compared a cooling flow solution directly to the observational data for the Phoenix cluster, rather than to a model fit of the data as done here. They found that a cooling flow provides a good match to the observations, supporting our conclusion that the ICM of this cluster may form a cooling flow.

5 Discussion

In the previous sections, we demonstrated that initially hydrostatic gaseous halos converge onto one of a single-parameter family of cooling flow solutions, within a cooling time. In these solutions the volume-filling gas phase flows inward at a velocity $v(r)\approx-r/t_{\rm cool}$, while the temperature is roughly equal to the circular temperature $T_{\rm c}\equiv\mu m_{\rm p}v_{\rm c}^{2}/(\gamma k_{\rm B})$. These solutions are similar to the solutions developed for the intracluster medium in the 1980’s (e.g. Fabian et al. 1984).

In Figure 14 we demonstrated that cooling flow solutions fail to explain the properties of the X-ray emission in typical clusters, as expected, though they are potentially consistent with the most dynamically-relaxed systems. On the other hand, we showed that a cooling flow solution with ${\dot{M}}\sim{\rm SFR}$ correctly predicts the X-ray absorption line properties of the Milky Way halo (Fig. 11), and the slope of the hot gas density profile derived by modelling the X-ray line emission in the MW halo (Fig. 12). Given that the cooling flow model is the simplest possible model and has no free parameters beyond the uncertainty in gas metallicity, this success may suggest that hot gas in the Milky Way halo forms a cooling flow with ${\dot{M}}\sim{\rm SFR}$, at least in the inner 10s of kpc which dominate the observed X-ray emission and absorption. However, Figure 13 shows that similar solutions underpredict the observed O vi by a factor of $\sim 5$, both in the MW halo and in the halos of external star-forming galaxies. The observations of O vi thus suggest that the cooling flow picture is incomplete.

In this section we compare cooling flow solutions with other existing models of gas in galaxy-scale halos, and discuss possible resolutions to the apparently inconsistent conclusions implied by the X-ray and O vi observations.

5.1 Comparison with hydrostatic models

Several authors have proposed hydrostatic models for gas in galaxy-scale dark matter halos (Maller & Bullock 2004; Faerman et al. 2017; Mathews & Prochaska 2017; Voit 2018). Since the hydrostatic equation on its own is insufficient to fully specify the structure of the halo gas, these studies evoked additional assumptions on the density or entropy profile. Specifically, Maller & Bullock (2004) assumed a flat entropy profile out to $R_{\rm cool}$, Faerman et al. (2017) derived the density profile from the average X-ray emission and absorption in the MW and from the O vi absorption around external galaxies, Mathews & Prochaska (2017) assumed a density which follows the dark matter density with a flat entropy core, and Voit (2018) assumed an entropy profile which yields $t_{\rm cool}/t_{\rm ff}\approx 10$ at all radii. In the cooling flow solutions for the MW halo gas (Figs. 11 – 13), the pressure profile is also close to hydrostatic, since the correction to the momentum equation is of order $\mathcal{M}^{2}$, while $\mathcal{M}$ is small at halo radii ($\approx 0.1$ at $100\,{\rm kpc}$, eqn. 32). However, although that the flow is relatively slow, it is this weak flow which sets the entropy profile of the halo gas via the entropy equation (3), at all radii smaller than $R_{\rm cool}$. Specifically, for Milky Way-mass halos we get $K\propto r^{0.87}$ (eqn. 25 with $m=-0.1$). Thus, in the absence of physical processes other than cooling, the entropy and density profiles in a quasi-static halo are not free parameters but rather set by the cooling flow solution.

5.2 Reconciling OVI and X-ray observations

Figure 15 summarizes the constraints on ${\dot{M}}$ derived in this study and in previous studies based on X-ray observations of O vii and O viii (black), and based on UV observations of O vi (blue). The horizontal position of each marker denotes the median distance contributing to the observations, while the range spans the $16-84$ percentiles. For the X-ray-based constraints this characteristic distance is estimated from the modelling, which for a density profile of $n_{\rm H}\propto r^{-1.6}$ predicted by cooling flows (eqn. 29) yields an X-ray gas distance of $\sim 10-40\,{\rm kpc}$. As discussed in section 4.1, comparison of the cooling flow solution with X-ray absorption yields ${\dot{M}}\approx 0.12-2.4$ times the MW SFR of $1.6\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ (black square in Fig. 15), while a comparison with the X-ray line emission yields a cooling rate of roughly half the SFR (black pentagon, see section 4.2 and Fig. 12). These deduced density profile, characteristic distance, and cooling rate of the X-ray gas are all consistent with the results of LB17 (black diamond), who modelled the X-ray observations using a phenomenological model independent of our cooling flow solutions. The results of LB17 are based on and consistent with other results by the same group (Miller & Bregman 2013, 2015; Bregman et al. 2018; Qu & Bregman 2018). The cooling flow solution hence provides a physical explanation for the parameters implied by LB17’s phenomenological analysis.

In contrast with the relative proximity implied for the X-ray gas, a deprojection of O vi observations around external galaxies shows that O vi absorption originates in gas farther out in the halo, at a radius range of $0.35-0.8R_{\rm vir}$, or $70-150\,{\rm kpc}$ for the median $R_{\rm vir}$ of $190\,{\rm kpc}$ in the COS-Halos sample (fig. 1 in S18, see consistent result in Mathews & Prochaska 2017). Cooling flow solutions with ${\dot{M}}={\rm SFR}$ underpredict the observed O vi absorption by a factor of $\sim 5$, and there is no stable cooling flow solution which can reproduce the observed O vi (Fig. 13, section 4.4). Fig. 15 plots previous estimates of the O vi-traced ${\dot{M}}$ assuming that it originates in radiatively-cooling gas, from Mathews & Prochaska (2017, blue star), from McQuinn & Werk (2018, blue cross) and from S18 (their ‘high-pressure’ scenario, blue plus). In all cases, we normalized the mass inflow rates found by these studies by the mean ${\rm SFR}=4.2\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ in the observed sample (Werk et al. 2013). The derived ${\dot{M}}/{\rm SFR}$ are roughly an order of magnitude larger than implied by the X-ray emission.

The O vi columns and dispersion observed around the MW are similar to those around external SF galaxies, if one accounts for the different perspectives and for confusion with the Galactic disk (Zheng et al. 2015, Fig. 13). It thus seems likely that O vi in the Milky Way has the same origin as O vi in other SF galaxies. Fig. 15 therefore suggests a somewhat surprising result that cooling rates in the O vi-bearing outer halo are significantly different from cooling rates in the inner halo which dominates the X-ray observations. It is important to note that this conclusion is already implied by the comparison of the LB17 and e.g. McQuinn & Werk (2018) data points even without considering the cooling flow analysis in this work, as discussed in Qu & Bregman (2018). Our cooling flow analysis strengthens this conclusion by providing a physical basis for LB17’s phenomenological results, as mentioned above. We next discuss several possible resolutions to this apparent tension between the X-ray and O vi observations.

Preferential heating of the OVI-bearing outer halo: To maintain the large cooling rates implied by the O vi observations, previous studies typically assumed that gas at $\gtrsim 100\,{\rm kpc}$ is heated by supernovae or AGN in the central galaxy (Cen 2013; Liang et al. 2016; Mathews & Prochaska 2017; Faerman et al. 2017; Suresh et al. 2017; McQuinn & Werk 2018). If this is the case, a heating mechanism powered by the central galaxy must avoid depositing most of its energy in the inner 10s of kpc to explain the order-of-magnitude lower cooling rates inferred from the X-ray observations. Alternatively, the heating energy may originate from larger radii, e.g. processes associated with cosmological structure assembly.

In the middle panel of Fig. 13 we plot a potential O vi profile in the presence of a heating mechanism that preferentially heats the outer halo (black dash-dotted line). To produce this profile we run a simulation similar to the simulations shown in Fig. 5, with $M_{\rm halo}=6\cdot 10^{11}\,{\rm M_{\odot}}$ and an initial baryon mass of $0.4$ times the halo baryon budget. At each assumed impact parameter, we integrate the O vi volume density in the $t=7\,{\rm Gyr}$ snapshot out to a maximum radius of $200\,{\rm kpc}$. This snapshot is chosen since it has a relatively large cooling rate of $20\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ within $200\,{\rm kpc}$, compared to a lower cooling rate of $4\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ within $40\,{\rm kpc}$, thus mimicking the different ${\dot{M}}$ implied by the O vi and X-ray observations shown in Fig. 15. Figure 13 demonstrates that in this snapshot the O vi columns are roughly consistent with observations. In our simulations which do not include heating mechanisms, such a gas structure is necessarily a transient phenomenon, which occurs only when ${\rm d}R_{\rm cool}/{\rm d}t$ increases quickly with time and there is hence a disparity between ${\dot{M}}$ at large scales and ${\dot{M}}$ at small scales (see eqn. 43). A transient solution is unlikely to explain the observed O vi which are seen around practically all blue galaxies. Potentially, in the presence of a heating mechanism which preferentially heats the outer halo as suggested here, such an O vi profile may be long-lived and thus consistent with the observations.

The low-pressure OVI scenario: The assumption that O vi traces radiatively cooling gas at $T\approx 10^{5.5}\,{\rm K}$ is mainly based on the results of cosmological simulations (Stinson et al. 2012; Hummels et al. 2013; Cen 2013; Oppenheimer et al. 2016; Liang et al. 2016; Gutcke et al. 2017; Suresh et al. 2017; Nelson et al. 2018). From a purely observational perspective, O vi along random quasar sight lines and in sight lines farther than $0.5R_{\rm vir}\approx 100\,{\rm kpc}$ from galaxies is almost always observed together with H i, with a column ratio of $N_{\rm HI}/N_{\text{\sc O\,vi}}\approx 1-10$. These H i observations support a scenario where the O vi and associated H i columns trace single-phase photoionized gas in thermal equilibrium with the UV background (Tripp et al. 2008; Thom & Chen 2008; Stern et al. 2016, S18). If the O vi-gas is in thermal equilibrium, then it is not radiatively cooling, and hence does not imply cooling rates significantly larger than suggested by the X-ray observations.

O vi can trace gas in thermal equilibrium if the gas has a relatively low thermal pressure of $\lesssim 1\,{\rm cm}^{-3}\,{\rm K}$, as expected outside the accretion shock. S18 argued that this would require a shock radius around \sim$$L^{\star} galaxies $R_{\rm shock}\lesssim 0.5R_{\rm vir}\approx 100\,{\rm kpc}$, substantially smaller than predicted by cosmological simulations. They showed that such a scenario is consistent with the observed C iii absorption and lack of low-ion absorption at $R_{\perp}\gtrsim 0.5R_{\rm vir}$, and with the observed linear relation between O vi column and velocity width. Also, if O vi traces free-falling preshock gas then the implied mass inflow rate is $\approx 5\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ (eqn. 30 in S18), comparable to the mean ${\rm SFR}=4.2\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ of the central galaxies. This alternative estimate of ${\dot{M}}/{\rm SFR}\approx 1.2$ implied by O vi is plotted in Fig. 15 (blue circle), and is comparable to the ${\dot{M}}/{\rm SFR}\sim 0.5$ suggested by the MW X-ray observations, thus alleviating the tension implied by assuming that O vi traces radiatively-cooling gas.

In the middle panel of Figure 13 we plot a possible $N_{\text{\sc O\,vi}}$ profile based on the model suggested by S18 (black dashed line), calculated as follows. For the hot gas inside the shock radius, we calculate a cooling flow solution with the mean $M_{\rm halo}=6\cdot 10^{11}\,{\rm M_{\odot}}$ and $z=0.2$ of the galaxies in the panel, ${\dot{M}}$ equal to their average SFR of $3\,{\rm M_{\odot}}\,{\rm yr}^{-1}$, and an assumed metallicity of $0.5\,{\rm Z_{\odot}}$. We require the outer boundary of the cooling flow to satisfy shock jump conditions with zero shock velocity. There is only a single possible solution for jump conditions consistent with preshock gas in thermal equilibrium with the UV background ($\approx 3\cdot 10^{4}\,{\rm K}$). This solution is plotted in Figure 18 in the appendix. For the gas outside the shock radius, the jump conditions imply the same ${\dot{M}}=3\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ as in the postshock gas and a preshock velocity222The characteristic velocity offset between the central galaxies and the O vi absorption profiles is $\approx 100\,\rm km\ s^{-1}$ (Tumlinson et al. 2011), lower than the deduced radial velocity of $150\,\rm km\ s^{-1}$ for the O vi gas in this model. This difference may potentially be explained by projection effects. Alternatively, assuming a preshock gas temperature somewhat higher than equilibrium with the background (e.g. due to adiabatic compression), would imply a lower preshock velocity which is more consistent with the observations. of $150\,\rm km\ s^{-1}$. We assume a velocity profile of $v\propto r^{-1/2}$ and a clumping factor of $\delta\rho/\rho=4$, which yield a best ‘by-eye’ fit to the O vi observations. This clumping factor is consistent with the O vi absorber pathlengths of $\gtrsim 10$s of kpc implied by observations relative to the total pathlength through the outer halo (Table 1 in S18). As can be seen in Figure 13, the expected $N_{\text{\sc O\,vi}}$ in this combined free-fall + cooling flow model is roughly consistent with the observed $N_{\text{\sc O\,vi}}$, while also matching the ${\dot{M}}\sim{\rm SFR}$ suggested by the MW X-ray observations. Predictions of this model for observations of Ne viii and other ions observable in the extreme UV (EUV) are discussed in Stern et al. (2018).

A density profile shallower than a cooling flow: An alternative resolution is that the halo density profile is shallower than deduced by LB17 from the MW observations (see section 2.4 in Bregman et al. 2018 for a discussion of this possibility). A flatter density profile would yield a larger characteristic radius and higher cooling rates for the X-ray gas. Such a flat density profile was deduced by Faerman et al. (2017), who fit a hydrostatic model to the average of the X-ray observations in the MW and to the O vi observations around external galaxies. The density and mass profile in their model are plotted in Figure 12. Note that Faerman et al. did not incorporate the dependence of X-ray observations on direction in the halo, which was used by LB17 to infer the slope of the density profile. In the Faerman et al. (2017) model, both the X-ray and O vi trace cooling rates of $\approx 30\,{\rm M_{\odot}}\,{\rm yr}^{-1}$. This model suggests both significant feedback heating in the halo and a closed baryon fraction (see lower panel of Figure 12), in contrast with the baryon-deficient halo and lack of heating implied by the cooling flow model and by LB17.

5.3 Comparison with precipitation models

We demonstrate in Figure 10 and section 3 that small amplitude perturbations in a cooling flow do not develop into multi-phase structure at radii where $t_{\rm cool}/t_{\rm ff}>1$, consistent with previous results based on linear perturbation theory (e.g. Balbus & Soker 1989) and hydrodynamic simulations (e.g. Joung et al. 2012; Sharma et al. 2012a). In this regime, the gas is predominantly single phase (by volume) and cool clouds develop only near strong disturbances which can seed non-linear perturbations (see section 4.5).

However, numerical simulations have demonstrated that in the presence of heating by feedback multi-phase structure may develop even if $t_{\rm cool}/t_{\rm ff}$ is as high as $\sim 10$ (Sharma et al. 2012a; Gaspari et al. 2013; Li & Bryan 2014b; Choudhury & Sharma 2016; Choudhury et al. 2019). Based on this result Sharma et al. (2012b) and Voit et al. (2015, 2017) proposed a feedback-regulated limit cycle. In this cycle, cool clumps in the halo ‘precipitate’ onto the central galaxy, providing fuel for feedback which heats the ambient medium and promotes further condensation. The main ansatz of these precipitation models, based on the results of numerical simulations, is that the feedback loop regulates the halo gas to have a minimum $t_{\rm cool}/t_{\rm ff}\sim 10$. Some support for this ansatz is provided by observations of X-ray emitting gas in groups and clusters of galaxies, which apparently adhere to the limit of $t_{\rm cool}/t_{\rm ff}\gtrsim 10$ and exhibit evidence for multiphase gas preferentially in objects and at radii where $t_{\rm cool}/t_{\rm ff}$ is close to the limit (McCourt et al. 2012; Voit & Donahue 2015).

In galaxy-scale halos, the roughly constant $t_{\rm cool}/t_{\rm ff}$ suggested by precipitation models is expected also in the absence of any feedback (eqn. 33 and lower-right panel of Figure 3). This result follows from the cooling flow solution (eqns. 24 – 27), which yields $t_{\rm cool}/t_{\rm ff}\propto r^{0.5+2m}$. For the weakly decreasing circular velocity profile expected in galaxy-scale halos ($v_{\rm c}\propto r^{m}\propto r^{-0.1}$, see Figure 1), we get $t_{\rm cool}/t_{\rm ff}$ which depends only weakly on radius. Thus, a roughly constant $t_{\rm cool}/t_{\rm ff}$ in simulations or implied by observations does not necessarily imply that galaxy-scale halos are precipitation-regulated.

In cooling flow solutions, the normalization of $t_{\rm cool}/t_{\rm ff}$ is set by the free parameter of the solution and in principle can have any value larger than unity. In practice, we expect the total hot gas mass – and thus the normalization of $t_{\rm cool}/t_{\rm ff}$ – to be set either by the cosmic baryon budget or modified by the loss of baryons driven by galactic outflows at high redshift. Interestingly, for a Milky Way-mass halo, assuming ${\dot{M}}=1\,{\rm M_{\odot}}\,{\rm yr}^{-1}\approx{\rm SFR}$ implies a characteristic $t_{\rm cool}/t_{\rm ff}=7.5$ (eqn. 33), approximately equal to the basic ansatz of precipitation models. Thus precipitation-regulated models for the MW yield roughly the same halo gas structure as the no-ongoing-feedback cooling flow solution discussed in this work. This similarity can be seen in the top panel of Figure 12, where the density profile of a precipitation-regulated halo (Voit 2018, black dash-dotted line) is compared to the density profile of the cooling flow solution (cyan line). The two profiles differ by less than a factor of two at all plotted radii, within the uncertainty of the modelling and of the observational constraints. Thus, the predictions of precipitation-regulated models for the hot gas structure and their comparison with observations (e.g. Voit 2018; Voit et al. 2019) do not differentiate between a precipitation-regulated CGM and a cooling flow CGM.

5.4 Additional implications for feedback

Figures 11–12 demonstrate that observational constraints on the hot gas in the MW halo are consistent with the predictions of a cooling flow with ${\dot{M}}\approx{\rm SFR}$. In this scenario ongoing heating by feedback is weak, and does not significantly change the halo gas properties. Strong feedback is however still required in this scenario at earlier epochs to explain the high metallicity ($0.3\,{\rm Z_{\odot}}-\,{\rm Z_{\odot}}$) and low gas mass ($\lesssim 20\%$ of the cosmic baryon budget, section 4.2) of the cooling flow solutions. Such a scenario where feedback is strong at high redshift but weak in the local Universe is qualitatively consistent with observational indications that galactic winds are strong in high SFR surface density galaxies common at $z\sim 2$, while they are weak in the low SFR surface density galaxies common at $z\sim 0$ (e.g., Heckman & Thompson 2017). This scenario is also (qualitatively) supported by cosmological simulations that predict strong winds at $z\sim 2$ which subside by $z\sim 0$ in $\sim L^{*}$ halos (e.g., in the FIRE cosmological simulations, Muratov et al. 2015, 2017).

Another requirement of the cooling flow model is that ongoing heating by the central black hole is inefficient. In the MW, possible evidence for such heating are the Fermi bubbles (Su et al. 2010), which extend $\approx 10\,{\rm kpc}$ from the Galactic center. Sarkar et al. (2017) estimated a feedback energy injection rate into the bubbles of $\approx(0.7-1)\cdot 10^{41}\,\rm erg\,{\rm s}^{-1}$. If this energy eventually heats the halo gas, it could offset the energy radiated by a cooling flow with ${\dot{M}}=1.5-2\,{\rm M_{\odot}}\,{\rm yr}^{-1}$ at radii $10-100\,{\rm kpc}$. Hence for feedback heating to be small relative to ${\dot{M}}\approx 1.6\,{\rm M_{\odot}}\,{\rm yr}^{-1}$, this heating rate estimate needs to be somewhat biased high, or alternatively the bubble energy is not efficiently transferred to the surrounding CGM.

In cluster-scale halos, the similarity of the entropy profiles of Phoenix and CLJ2043-5035 with cooling flow solutions may suggest that these clusters are in a cooling flow phase of a feedback limit cycle (see introduction), as also supported by the large UV-estimated SFRs in their central galaxy ($\approx 2000$ and $\approx 200\,{\rm M_{\odot}}\,{\rm yr}^{-1}$, respectively, McDonald et al. 2012, 2016). It would be useful to conduct a similar comparison of the cooling flow solutions with simulations of feedback limit cycles (e.g. Prasad et al. 2015) and thus derive the predicted duty cycle of the cooling flow phase in such models. This prediction could then be compared to the observed fraction of clusters in the cooling flow phase.

6 Summary

In this paper, we solve the spherical steady-state equations of radiatively cooling gas, in a time-independent gravitational potential characteristic of dark matter halos. Our solutions extend previous solutions derived in the 1980’s for cluster-scale halos to galaxy-scale halos. We derive self-similar solutions in the hydrostatic limit for power-law potential and cooling function, and numerical solutions for more general potentials and cooling functions. Both type of solutions have a single free parameter, and are thus fully determined if either the mass inflow rate or total halo gas mass are known. We find that in cooling flows the entropy scales with radius as $K\propto r^{1+4m/3}$ (eqn. 13), where $m\equiv{\rm d}\log v_{\rm c}/{\rm d}\log r$ is the power-law index of the circular velocity radial profile.

Using idealized 3D hydrodynamic simulations of initially hydrostatic gas, we showed that the average gas properties in the simulation converge onto a cooling flow solution within a cooling time or, equivalently, out to the cooling radius (Figs. 6 – 9). Density fluctuations remain significantly below unity as long as $t_{\rm cool}>t_{\rm ff}$ (Fig. 10), consistent with expectations based on linear theory. Halo gas which forms a cooling flow is thus expected to be predominantly single phase, with cool clouds limited to strong disturbances that can seed non-linear perturbations (e.g. near the disk plane or due to stirring by satellite galaxies).

We compare the cooling flow solutions with observational constraints on halo gas at low redshift. Our conclusions can be summarized as follows:

Observations of O vii and O viii absorption in the Milky Way halo are consistent with a cooling flow solution with ${\dot{M}}\sim{\rm SFR}$, with no free parameters (Fig. 11). 2. 2.

The Milky Way halo gas density profile of $n_{\rm e}\propto r^{-1.5}$, deduced by Li & Bregman (2017) from modelling of observed O vii and O viii emission, is consistent with the prediction of cooling flows (Fig. 12). 3. 3.

Estimates of the thermal pressure based on cool clouds in the ISM and in the Magellanic Stream, and the dispersion measure towards the LMC, are also consistent with a cooling flow solution with ${\dot{M}}\sim{\rm SFR}$. 4. 4.

Cooling flows with ${\dot{M}}\sim{\rm SFR}$ underpredict observed O vi absorption columns in the MW halo and in the halos of other star-forming galaxies, typically by a factor of $\sim 5$, if the cooling flows are assumed to extend out to $\gtrsim R_{\rm vir}$ (Fig. 13). This discrepancy can be reconciled with the successes of the cooling flow solution described in (i)–(iii) if the O vi-bearing outer halo is preferentially heated or, alternatively, if the large $N_{\text{\sc O\,vi}}$ originate in free-falling gas outside an accretion shock at $\lesssim 100\,{\rm kpc}$ which subsequently develops into a cooling flow at smaller radii, as suggested by Stern et al. (2018). Upcoming observations of Ne viii and other ions observable in the EUV could potentially test these possibilities. 5. 5.

The cooling flow solution for the MW halo gas predicts $t_{\rm cool}/t_{\rm ff}\approx 7.5(r/100\,{\rm kpc})^{0.3}$, similar to the basic ansatz of thermal-instability-regulated feedback-loop models, in which $t_{\rm cool}/t_{\rm ff}\approx 10$ independent of radius. These feedback-loop models thus yield roughly the same hot gas structure and observables as the no-feedback cooling flow solution (Fig. 12). 6. 6.

In galaxy cluster halos, cooling flows predict a radial entropy profile of $K\propto r^{1.4}$, steeper than the observed profile of $K\propto r^{0.7}$ in typical cool-core clusters. This is consistent with other evidence that feedback alters the gas properties in the intracluster medium. However, two of the most dynamically-relaxed clusters in the SPT sample have entropy profiles consistent with cooling flows beyond the galaxy scale $\approx 30\,{\rm kpc}$ (Fig. 14). These clusters may be in the cooling flow phase of a feedback limit cycle.

The possibility that hot gas in Milky-Way-like halos forms a cooling flow with negligible ongoing feedback heating, could be a consequence of the fact that galactic winds are weak in galaxies with low SFR surface density (e.g., Heckman & Thompson 2017). This is also suggested by simulations that predict weak winds in low-redshift \sim$$L^{\star} galaxies (e.g., Muratov et al. 2015), as discussed in section 5.4. Our analysis thus raises the question of which halos (as a function of mass and redshift) can be adequately modeled by a cooling flow over an interesting range of radii. At high redshift, stellar feedback is expected to be stronger and is also required to explain the high metallicity and low gas mass implied by the cooling flow solutions at low redshift. Comparison of cooling flow solutions with observations of high-redshift halos (e.g. Turner et al. 2014; Rudie et al. 2019) may be able to address this important question.

Acknowledgements

We thank the referee, Prateek Sharma, for a thorough report and illuminating comments that significantly improved the paper. We thank also Mark Voit, Sean D. Johnson, and Zachary Hafen for detailed and insightful comments. We thank Michael McDonald for the entropy profiles fits used in Fig. 14. JS is supported by the CIERA Postdoctoral Fellowship Program. DF is supported by the Flatiron Institute, which is supported by the Simons Foundation. CAFG is supported by NSF through grants AST-1517491, AST-1715216, and CAREER award AST-1652522, by NASA through grants NNX15AB22G and 17-ATP17-0067, by STScI through grants HST-GO-14681.011, HST-GO-14268.022-A, and HST-AR-14293.001-A, and by a Cottrell Scholar Award from the Research Corporation for Science Advancement. This work was supported in part by a Simons Investigator Award from the Simons Foundation and by NSF grant AST-1715070.

Appendix A Integration details

In this section we describe the numerical integration technique used to derive the marginally-bound transonic solutions.

We first integrate the flow equations outward from a sonic point at an assumed radius $R_{\rm sonic}$, in order to derive the subsonic part of the solution. In the integration we use the form of the flow equations expressed in eqns. (1), (13), and (17). The temperature at $R_{\rm sonic}$ is initially assumed to equal $T_{\rm c}$, while $\rho(R_{\rm sonic})$ and $v(R_{\rm sonic})$ are derived from this assumption and from the requirement that both sides of eqn. (17) vanish. Since the momentum equation is indeterminate at $R_{\rm sonic}$, we start the integration at a radius $R_{0}=(1+\epsilon)R_{\rm sonic}$ with $\epsilon=10^{-5}$. The values of the hydrodynamic variables $T(R_{0})$, $v(R_{0})$, and $\rho(R_{0})$ thus differ from their values at $R_{\rm sonic}$ by factors of $1+({\rm d}\ln T/{\rm d}\ln r)\epsilon$, $1+({\rm d}\ln v/{\rm d}\ln r)\epsilon$, and $1-({\rm d}\ln v/{\rm d}\ln r+2)\epsilon$, respectively. The derivation of ${\rm d}\ln T/{\rm d}\ln r$ and ${\rm d}\ln v/{\rm d}\ln r$ for a given $T(R_{\rm sonic})$ is described below. If the integrated solution reaches $\mathcal{M}=1$ with an infinite velocity derivative, then $T(R_{\rm sonic})$ is adjusted upward, and the flow equations are integrated again. If the solution becomes unbound ($\mathcal{B}>0$), then $T(R_{\rm sonic})$ is adjusted downward and again the integration is repeated. The process continues until we find a solution which is bound out to a large radius of $10R_{\rm vir}$. Then, using $T(R_{\rm sonic})$ of this marginally-bound solution we find also the supersonic part of the solution, by applying a similar offset to the hydrodynamic variables with $\epsilon=-10^{-5}$, and integrating inward. An example of this process is shown in Figure 16. The derived marginally-bound solution in this example is the solution shown in Fig. 2.

To derive ${\rm d}\ln T/{\rm d}\ln r$ and ${\rm d}\ln v/{\rm d}\ln r$ at the sonic point we divide eqn. (17) by $\mathcal{M}^{2}-1$, which gives

[TABLE]

For the velocity derivative to be finite at the sonic radius (where $v=-c_{\rm s}$), the numerator must equal zero, i.e.

[TABLE]

where for convenience we defined the parameter $x\equiv v_{\rm c}^{2}/2c_{\rm s}^{2}=T_{\rm c}/2T$. Note that for adiabatic conditions ($t_{\rm flow}/t_{\rm cool}\rightarrow 0$), condition (47) is equal to the standard Bondi condition of $v_{\rm c}^{2}=2c_{\rm s}^{2}$ or $x=1$. The logarithmic derivative of $v$ at the sonic point is then found using l’Hospital’s rule:

[TABLE]

The denominator is equal to

[TABLE]

where in the second equality we used $-v=c_{\rm s}$. The relation between the logarithmic derivative of $v$ and the logarithmic derivative of $T$ can be derived from mass and entropy conservation (eqns. 12 and 13):

[TABLE]

Using the sonic point equality (47) for $t_{\rm flow}/t_{\rm cool}$, assuming $\gamma=5/3$ and rearranging we get

[TABLE]

Plugging this equation in eqn. (49) then yields for the denominator of eqn. (48)

[TABLE]

For the numerator in eqn. (48) we note that the term $c_{\rm s}^{2}t_{\rm flow}/\gamma t_{\rm cool}$ is proportional to $r\rho\Lambda/-v$. Hence for $v_{\rm c}=v_{\rm c}(r)$ and $\Lambda=\Lambda(T,\rho)$ we get

[TABLE]

Using again mass conservation (${\rm d}\ln\rho/{\rm d}\ln r=-{\rm d}\ln v/{\rm d}\ln r-2$) and the sonic point equality (47) we get

[TABLE]

Using eqn. (51) in eqn. (A) and rearranging we get for the numerator of eqn. (48)

[TABLE]

Finally, using equations (51), (52), and (A) in eqn. (48) we get that

[TABLE]

which can be solved as a quadratic equation for a given $x=v_{\rm c}^{2}/2c_{\rm s}^{2}$, ${\rm d}\ln\Lambda/{\rm d}\ln T$, ${\rm d}\ln\Lambda/{\rm d}\ln\rho$, and ${\rm d}\ln v_{\rm c}/{\rm d}\ln r$. Typically, one of these roots has ${\rm d}\mathcal{M}/{\rm d}r<0$, which corresponds to a cooling flow solution that is subsonic at large scales and supersonic at small scales, while the other root has ${\rm d}\mathcal{M}/{\rm d}r>0$, which corresponds to the opposite transition. We use the former root in the solutions used in this paper.

A similar derivation was done by Mathews & Guo (2012) for the specific case of a point mass and constant cooling function, and their result can be reproduced using ${\rm d}\ln\Lambda/{\rm d}\ln T={\rm d}\ln\Lambda/{\rm d}\ln\rho=0$ and ${\rm d}\ln v_{\rm c}/{\rm d}\ln r=-1/2$ in eqn. (56). Note that the parameter ‘a’ defined by Mathews & Guo is equal to $(1-x)/x$ in our notation.

Appendix B Outer boundary conditions

The outer boundary condition of the transonic marginally-bound solutions used in this work is $\mathcal{B}\rightarrow 0^{-}$ as $r\rightarrow\infty$ (see Appendix A). How will the solution change if we impose a different outer boundary condition? In Figure 17 we repeat the integration process shown in Fig. 16, but applying the shooting method until some specific temperature at $200\,{\rm kpc}$ is reached. The marginally-bound solution is marked by a thick black line. For the alternative boundary conditions, we choose only temperatures which yield a bound ($\mathcal{B}<0$) outer boundary condition, since gas with $\mathcal{B}>0$ will likely escape from the halo rather than form an accretion flow. As can be seen in the plot, the solutions deviate from the marginally-bound solution only near the outer boundary, with a maximum deviation in temperature at $100\,{\rm kpc}$ of $\pm 0.15\,{\rm dex}$. A similar conclusion is reached if we impose some gas density at the outer boundary. The exact choice of the outer boundary condition thus does not affect our conclusions.

How would the solution change if the outer boundary is an accretion shock? In Figure 18 we compare a marginally-bound solution with the solution discussed in section 5.2, which satisfies shock jump conditions at $R_{\rm shock}=100\,{\rm kpc}$. Both solutions have the same ${\dot{M}}$. The solutions differ mainly at $\lesssim 100\,{\rm kpc}$, consistent with the conclusion from Fig. 17.

Appendix C The timescale for establishing the $t_{\rm cool}\approx t_{\rm flow}$ relation

Figures 6 – 9 show that in the hydrodynamic simulations $t_{\rm cool}\approx t_{\rm flow}$ out to radii $r\gg R_{\rm cool}$, indicating that this relation is established on a timescale $\ll t_{\rm cool}$. In this section we explain the origin of this behavior.

At times $t\ll t_{\rm cool}$, the radial momentum equation of an initially hydrostatic halo is

[TABLE]

where $P_{0}$ is the gas pressure at $t=0$ which satisfies $\rho^{-1}\partial P_{0}/\partial r=v_{\rm c}^{2}/r$. The flow acceleration hence equals

[TABLE]

where in the last ‘$\sim$’ we used $r/c_{\rm s}\sim t_{\rm ff}$ and disregarded order-unity factors. At $t\approx t_{\rm ff}$, the term on the right is of order $r/(t_{\rm cool}t_{\rm ff})$, indicating that a velocity of order $r/t_{\rm cool}$ can be reached within a free-fall timescale. The $t_{\rm cool}\approx t_{\rm flow}$ relation is thus established on a sound-crossing timescale rather than on a cooling timescale.