A Global Model for Circumgalactic and Cluster-core Precipitation

All modern discussions of thermal instability in the context of galaxy formation are rooted in the classic works of Rees & Ostriker (1977), Binney (1977), and Silk (1977), which themselves owe a debt to Hoyle (1953). These landmark papers derived the maximum stellar mass of a galaxy ($\sim {10}^{12}\,{M}_{\odot }$) by comparing the time for gas to fall through a galaxy's potential well to the time required to cool from the potential's virial temperature. If the cooling time is less than the freefall time, then infalling gas can potentially condense and fragment into star-forming clouds via thermal instability. That can happen relatively easily in galaxy-scale objects with virial temperatures $\lesssim {10}^{7}$ K but is more difficult to arrange in hotter, more massive systems, leading to a natural division between the mass scales of individual galaxies and those of galaxy groups and clusters.

Many subsequent papers have made interesting use of the cooling time to freefall time ratio to analyze how galaxies form and evolve (e.g., Blumenthal et al. 1984; Fall & Rees 1985; Maller & Bullock 2004). However, some form of negative feedback is necessary to explain the inefficient transformation of a galaxy's gas supply into stars (e.g., Larson 1974; White & Rees 1978; Dekel & Silk 1986; White & Frenk 1991; Baugh et al. 1998; Somerville & Primack 1999; Kauffmann & Haehnelt 2000). A complete understanding of the relationship between thermal instability and galaxy formation must therefore account for the interplay between radiative cooling and the energetic feedback that opposes cooling.

Field (1965) worked out many of the fundamental features of astrophysical thermal instability but did not consider the complications that arise when thermal instability couples with buoyancy. The most comprehensive analytical treatment of that coupling is by Balbus & Soker (1989), who astutely summarized much of the preceding work. Balbus & Soker were primarily concerned with the development of inhomogeneity in galaxy-cluster cores. At the time, hot gas in the cores of many galaxy clusters was suspected to condense at rates $\sim {10}^{2-3}\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ (e.g., Fabian 1994), but models of homogeneous cooling flows into the cluster's central galaxy produced X-ray surface-brightness profiles with central peaks far greater than what were observed. This conundrum led to speculation that the mass inflow rate could decline inward because of thermal instability and spatially distributed condensation, which would reduce the radiative losses required to maintain a steady state at the center of the flow (e.g., Thomas et al. 1987; White & Sarazin 1987a, 1987b, 1987c). Upon closer examination, this speculation was found to be problematic, because buoyancy generally tends to suppress the development of thermal instability (e.g., Cowie et al. 1980; Nulsen 1986; Balbus & Soker 1989). The attention of the field therefore gradually shifted away from steady-state cooling-flow models in favor of models in which feedback from a central active galactic nucleus compensates for cooling (e.g., Tabor & Binney 1993; Binney & Tabor 1995; Soker et al. 2001; McNamara & Nulsen 2007, 2012).

More than two decades later, the coupling between buoyancy and thermal instability is being reexamined, because the presence of inhomogenous gas in galaxy-cluster cores consisting of multiple phases that are orders of magnitude cooler and denser than the ambient medium now appears closely linked with the ratio of cooling time to freefall time (e.g., Gaspari et al. 2012; McCourt et al. 2012; Voit & Donahue 2015; Voit et al. 2015b). In those studies, the cooling time is typically defined with respect to the specific heat at constant volume, so that ${t}_{\mathrm{cool}}\equiv [3{kT}/{n}_{e}{\rm{\Lambda }}(T)](n/2{n}_{i})$, where ${\rm{\Lambda }}(T)$ is the usual cooling function at temperature T, and the number densities of electrons, ions, and gas particles are n_e , n_i , and n, respectively. The freefall time ${t}_{\mathrm{ff}}\equiv {(2r/g)}^{1/2}$ is defined with respect to the local gravitational potential g at radius r. Given these definitions, a floor appears to be present near ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\approx 10$ in the radial cooling-time profiles of the ambient hot gas in galaxy clusters. A large majority of the cluster cores known to contain multiphase gas have minimum values of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ within a factor of 2 of this floor. Conversely, almost all of the clusters without multiphase gas have $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\gt 20$.

McCourt et al. (2012) and Sharma et al. (2012b) interpreted this relationship between multiphase gas and ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ as resulting from the amplification of initially small perturbations by thermal instability. Under conditions of global thermal balance, numerical simulations of thermal instability in a plane-parallel potential by McCourt et al. (2012) showed that ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\lesssim 1$ was required for thermal instability resulting in condensation, but the simulations of Sharma et al. (2012b) in a spherically symmetric potential indicated that condensation could happen in spherical systems with ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\lesssim 10$. The threshold value of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ was therefore assumed to be geometry dependent. Subsequent simulations implementing more sophisticated treatments of feedback appeared to corroborate that assumption because they showed that condensation-fueled accretion of cold gas onto a central black hole can lead to long-lasting self-regulation with $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$ (e.g., Gaspari et al. 2012; Li & Bryan 2014a; Li et al. 2015; Prasad et al. 2015).

In the meantime, it has become clear that the critical ratio ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ is not a geometry-dependent manifestation of local thermal instability. For example, numerical simulations by Meece et al. (2015) of thermal instability in a plane-parallel potential under conditions seemingly quite similar to those adopted by McCourt et al. (2012) showed that condensation could occur at the midplane of systems with any value of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$. Also, Choudhury & Sharma (2016) presented a detailed thermal stability analysis of systems in global thermal balance showing that the growth rates of linear perturbations are largely independent of the gravitational potential's geometry. So why, then, do both real and simulated galaxy-cluster cores appear to self-regulate at $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$?

The most general answer seems to involve a phenomenon that we have come to call precipitation. As feedback acts on a galaxy-scale system, the outflows it drives can promote condensation of the hot ambient medium by raising some of it to greater altitudes (e.g., Revaz et al. 2008; Li & Bryan 2014b; McNamara et al. 2016). Adiabatic uplift promotes condensation by lowering the ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ ratio of the uplifted gas (see Section 2). This process is loosely analogous to the production of raindrops during adiabatic cooling of uplifted humid gas in a thunderstorm—hence, the name "precipitation." As in a thunderstorm, the condensates rain down toward the bottom of the potential well after they form. In simulations, this rain of cold gas into the galaxy at first provides additional fuel for feedback and temporarily boosts the strength of the outflows, but eventually those strengthening outflows add enough heat to the ambient medium to raise ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ high enough to stop the condensation. Precipitation is therefore naturally self-regulating.

A prescient series of "cold feedback" papers by Soker & Pizzolato anticipated many features of the precipitation mechanism now seen in simulations (Pizzolato & Soker 2005, 2010; Soker 2006, 2008). They proposed a feedback cycle in cluster cores in which energetic AGN outbursts produce a wealth of nonlinear density perturbations, the densest of which cool, condense, fall back toward the black hole, and provide more fuel for accretion. In this scenario, the cooling times of the blobs must be short enough that ${t}_{\mathrm{cool}}\lesssim {t}_{\mathrm{ff}}$ and also shorter than the time interval between AGN heating outbursts. They argued that this source of accretion fuel could potentially provide much more fuel than Bondi accretion from the hot medium alone while responding to changes in the ambient cooling time far more quickly. They also noted that a shallow central entropy gradient would promote condensation.

Shortly thereafter, Gaspari and collaborators produced sets of numerical simulations in which such a cold feedback loop was realized (Gaspari et al. 2012, 2013, 2015; but see also Sharma et al. 2012b). In these simulations, cold clouds form through thermal instability and accrete toward the center of the simulation volume while the heating mechanism needed to maintain approximate global thermal balance stimulates turbulence. This turbulence is critical, because it ensures a steady supply of cold gas blobs with low specific angular momentum, which can plunge to the center through a process the authors call "chaotic cold accretion."

The proliferation of terminology has a way of making these mechanisms seem more different than they really are. Precipitation-regulated feedback, as described in this paper, is a "cold feedback" mechanism that fuels a central black hole through "chaotic cold accretion." Here we are proceeding with the term "precipitation," despite the prior existence of these other terms, because the processes that promote thermal instability and condensation in circumgalactic gas may be responsible for more than just the feeding of black holes. Therefore, they merit a more general term.

Precipitation is potentially of broader interest because it may link a galaxy's time-averaged star formation rate with the multiphase structure of its circumgalactic gas. Observations of circumgalactic absorption lines are showing that the masses of gas and metals within a few hundred kiloparsecs of a galaxy are at least as great as the galaxy's mass in the form of stars (e.g., Tumlinson et al. 2011). Also, the amount of circumgalactic gas at intermediate temperatures (10⁵–10⁶ K) appears closely linked with a galaxy's star formation rate. Observations of the circumgalactic medium at other wavelengths likewise show rich multiphase structure (e.g., Putman et al. 2012), which has been challenging for simulations to reproduce (e.g., Hummels et al. 2013; Ford et al. 2016).

Sharma et al. (2012a) proposed that thermal instability, through precipitation-regulated feedback, places a lower limit of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\approx 10$ on the ambient density of circumgalactic gas (but see also Meece et al. 2015). Voit et al. (2015a) built on that idea to show how precipitation-regulated feedback could be responsible for governing not only the relationships between stellar mass, metallicity, and the stellar baryon fraction observed among galaxies but also the relationship between a galaxy's stellar velocity dispersion and the mass of its central black hole. However, they did so without having a satisfactory explanation for the crucial assumption that $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$ or a complete model for the global structure of precipitation-regulated systems.

This paper's purpose is to propose a global context for interpreting observations of multiphase gas around galaxies, as well as the numerical simulations that strive to reproduce those observations, in terms of the long legacy of theoretical papers on astrophysical thermal instability. Many of the theoretical results derived here were published decades ago by others, most notably by Defouw (1970), Cowie et al. (1980), Nulsen (1986), Malagoli et al. (1987), Loewenstein (1989); and Balbus & Soker (1989). Our rederivations of them are intended to provide a common conceptual framework for a meta-analysis of the simulations in which precipitation occurs. We are focusing on simulations of precipitation in galaxy-cluster cores, because that is where the observations of multiphase gas around galaxies and interactions of outflows with the circumgalactic medium are most complete. However, many of the physical principles that govern condensation in those environments apply to circumgalactic gas around galaxies of all masses.

Busy readers may wish to be selective in deciding which sections of this long paper will reward their close attention. For them, we have prepared this guide, along with a cartoon (Figure 1) that sketches out the main ideas.

Zoom In Zoom Out Reset image size

• 1.  
  The next three sections (Sections 2–4) consider the local conditions required for condensation and precipitation. Two different modes of precipitation emerge from those considerations. One is analogous to rain that is stimulated by the uplift of humid gas in Earth's atmosphere, because of the role that adiabatic cooling plays in bringing on condensation. The other is analogous to drizzle or fog, in that condensation happens without uplift, when the conditions are right.

  • (a)  
    Section 2 initiates the discussion with a brief reminder about the deep connections between adiabatic uplift and condensation. It also presents a short calculation showing that ambient gas uplifted at speeds comparable to a halo's circular velocity is likely to condense if it initially has ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\lesssim 10$. This finding suggests that the ambient medium around a galaxy cannot persist in a state with ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\ll 10$ if there is significant vertical circulation. Galactic outflows in which the energy source is fueled by condensation therefore tend to drive the ambient medium toward ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\approx 10$.
  • (b)  
    Section 3 outlines the general conditions necessary for thermal instability without uplift to progress to condensation in a hydrostatic medium that is thermally balanced within each equipotential layer. The main result is that the condition for condensation to occur depends not only on ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ but also on the slope of the entropy gradient: thermal instability leads to condensation only if a low-entropy perturbation can cool faster than it sinks to a layer of equivalent entropy. This is not a new result, but its significance is often not fully appreciated.
  • (c)  
    Section 4 uses the results of Section 3 to interpret recent simulations of thermal instability in circumgalactic gas. In particular, it calls attention to the critical role of the global entropy gradient in determining where condensation can occur and where it is suppressed by buoyancy. The main result is that media in thermal balance are prone to condensation in regions where the large-scale entropy gradient is flat. However, buoyancy tends to delay the onset of condensation if ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\gg 1$.
• 2.  
  The following two sections (Sections 5 and 6) apply the findings from the first part of the paper to interpret the global evolution of simulated galactic systems in which condensation fuels feedback. Our objective is to understand why those systems end up in long-lasting configurations with $5\lesssim \min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\lesssim 20$ and with radial entropy profiles in agreement with observations of multiphase galaxy-cluster cores.

  • (a)  
    Section 5 examines recent simulations of condensation in globally balanced but locally unstable galactic systems in light of the findings summarized in Section 4. It points out that feedback is required for the development of a multiphase medium, because phase separation cannot happen in a globally balanced medium without the flow of free energy through the system. Condensation also cannot happen through linear thermal instability in regions with a significant entropy gradient and ${t}_{\mathrm{cool}}\gg {t}_{\mathrm{ff}}$. Steady self-regulation of a precipitating system therefore favors a global configuration with a shallow inner entropy gradient and a steeper outer entropy gradient (Figure 1). The shallow inner gradient promotes the precipitation needed for fuel, while the steeper outer entropy gradient prevents condensation from running away into a cooling catastrophe. The boundary between these regions tends to be where ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ reaches a minimum value ∼10, for reasons outlined in Section 2. In order to ensure long-term global stability with a central cooling time $\lesssim 1\,\mathrm{Gyr}$, much of the feedback energy must propagate beyond the isentropic zone before thermalizing. This finding has deep implications for implementations of black-hole feedback in numerical simulations.
  • (b)  
    Section 6 puts all of these pieces together to interpret the time-dependent behavior of precipitation-regulated feedback in a numerical simulation from Li et al. (2015). It shows that unopposed cooling leads to a power-law entropy gradient at the center of the system, which focuses condensation onto the central black hole. The feedback response disrupts that central gradient out to where ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\approx 10$ in the ambient medium. Much of the gas uplifted from that central region can be induced to condense, particularly if it is inhomogeneous. The system then settles into a long-lasting steady state in which condensed gas fuels the outflow and feedback maintains the isentropic central region at a level corresponding to $5\lesssim \min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\lesssim 20$. Catastrophic cooling is prevented because much of the outflow's energy is thermalized outside of the isentropic zone. When the condensed gas is depleted, the outflow shuts down. Cooling then proceeds almost homogeneously, because there is no source of free energy to promote phase separation. However, thermal instability eventually initiates condensation near the outer edge of the isentropic region, at the minimum of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ in the ambient medium. Newly condensed gas subsequently reignites feedback, and the cycle repeats.

The paper concludes with two sections that acknowledge the loose ends (Section 7) and present some concluding thoughts about how the overall model can be tested (Section 8), followed by an Appendix that constructs a useful toy model for global configuration changes of precipitation-regulated systems.

To illustrate this point, consider a blob of gas that begins at radius r₁ with entropy K₁ in an isothermal potential with constant circular velocity v_c . The freefall time in this potential is linearly proportional to radius, and the cooling time of the blob depends on its pressure. When written in terms of the entropy K and gas pressure P, this cooling time is

$\begin{eqnarray}&&{t}_{\mathrm{cool}}=\displaystyle \frac{3{K}^{6/5}}{{P}^{1/5}\,{\rm{\Lambda }}[T(P,K)]}{\left(\displaystyle \frac{n}{{n}_{e}}\right)}^{1/5}\left(\displaystyle \frac{n}{2{n}_{i}}\right)\end{eqnarray} \tag{ 1 }$

$\begin{eqnarray} & \approx & 235\,\mathrm{Myr}\,{\left(\displaystyle \frac{K}{10\mathrm{keV}{\mathrm{cm}}^{2}}\right)}^{6/5}{\left(\displaystyle \frac{P}{k\cdot {10}^{6}{\rm{K}}{\mathrm{cm}}^{-3}}\right)}^{-1/5}\\ & & \times {\left(\displaystyle \frac{{\rm{\Lambda }}}{2\times {10}^{-23}\mathrm{erg}{\mathrm{cm}}^{3}{{\rm{s}}}^{-1}}\right)}^{-1}.\end{eqnarray} \tag{ 2 }$

For gas with the temperatures and abundances typical of galaxy-cluster cores and large elliptical galaxies (${10}^{7}\,{\rm{K}}\lesssim T\lesssim 3\times {10}^{7}\,{\rm{K}}$), the cooling function ${\rm{\Lambda }}(T)$ is nearly independent of temperature, meaning that ${t}_{\mathrm{cool}}$ changes very little as gas adiabatically expands or contracts. One can account for the temperature dependence of Λ by defining $\lambda \equiv d\mathrm{ln}{\rm{\Lambda }}/d\mathrm{ln}T$ and noting that $T\propto {K}^{3/5}{P}^{2/5}$, to obtain ${t}_{\mathrm{cool}}\,\propto {K}^{(6-3\lambda )/5}{P}^{-(1+2\lambda )/5}$. Adiabatic uplift of a gas blob in a stratified atmosphere therefore changes the ratio of cooling time to freefall time within the blob according to

$\begin{eqnarray}&&{\left(\displaystyle \frac{{t}_{\mathrm{cool}}}{{t}_{\mathrm{ff}}}\right)}_{r}={\left(\displaystyle \frac{{t}_{\mathrm{cool}}}{{t}_{\mathrm{ff}}}\right)}_{{r}_{1}}{\left(\displaystyle \frac{{r}_{1}}{r}\right)}^{1-2\beta (1+2\lambda )/5},\end{eqnarray} \tag{ 3 }$

where $d\mathrm{ln}P/d\mathrm{ln}r=-2\beta$ with $\beta \equiv \mu {m}_{p}{v}_{c}^{2}/2{kT}$ in a hydrostatic ambient medium. Uplift by an order of magnitude in radius consequently lowers the ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ ratio local to the blob by roughly an order of magnitude, given the typical β and λ values of circumgalactic gas ($0.5\lesssim \beta \lesssim 1$ and $-1\lesssim \lambda \,\lesssim 0.5$), as long as magnetic fields sufficiently suppress thermal conduction.

This relationship helps explain why precipitation-driven feedback tends to push a galaxy's volume-filling ambient medium into a state with $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$. If the central gas is so dense that ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\ll 10$, then outflows propelled by feedback can promote additional condensation of the ambient medium with relatively modest amounts of adiabatic uplift. Newly condensed clouds that have not achieved escape velocity can then rain back down into the galaxy to fuel more star formation and AGN feedback, and the rain will continue until the dissipation of feedback energy raises the ambient value of $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})$ enough to inhibit condensation of uplifted gas. Then, feedback finally succeeds in stopping the accumulation of cold clouds within the galaxy. The critical value of $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})$ at which cessation of condensation should happen is hard to calculate from first principles, but simulations of the process indicate that it happens around $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$ (e.g., Gaspari et al. 2012; Li et al. 2015; Prasad et al. 2015; Meece et al. 2016).

A more quantitative relationship between uplift and condensation can be derived from the cooling history of an uplifted gas blob that follows a purely radial ballistic trajectory without incurring any heat input. For simplicity, consider a gravitational potential in which v_c is constant with radius, so that $\phi (r)={v}_{c}^{2}\mathrm{ln}(r/{r}_{\max })$, where ${r}_{\max }$ is the radius at which the blob reaches its apex. This implies

$\begin{eqnarray}&&r={r}_{\max }\exp \left(-\displaystyle \frac{{v}^{2}}{2{v}_{c}^{2}}\right),\end{eqnarray} \tag{ 4 }$

where v is the radial velocity, which leads to the equation of motion

$\begin{eqnarray}&&\ddot{r}=\dot{v}=-\displaystyle \frac{{v}_{c}^{2}}{{r}_{\max }}\exp \left(\displaystyle \frac{{v}^{2}}{2{v}_{c}^{2}}\right).\end{eqnarray} \tag{ 5 }$

We are interested in knowing the blob's time of flight ${\rm{\Delta }}t$ from an initial velocity v₁ to a final velocity v₂, which can be determined by integrating the equation of motion:

$\begin{eqnarray}{\rm{\Delta }}t & = & \displaystyle \frac{{r}_{\max }}{{v}_{c}^{2}}{\displaystyle \int }_{{v}_{2}}^{{v}_{1}}\exp \left(\displaystyle \frac{{v}^{2}}{2{v}_{c}^{2}}\right){dv}\\ & = & \displaystyle \frac{{\pi }^{1/2}}{2}\left[\mathrm{erf}\left(\displaystyle \frac{{v}_{1}}{\sqrt{2}{v}_{c}}\right)-\mathrm{erf}\left(\displaystyle \frac{{v}_{2}}{\sqrt{2}{v}_{c}}\right)\right]{t}_{\mathrm{ff}}({r}_{\max }).\end{eqnarray} \tag{ 6 }$

The maximum time of flight is then $\max ({\rm{\Delta }}t)={\pi }^{1/2}{t}_{\mathrm{ff}}({r}_{\max })$, and setting ${v}_{1}=-{v}_{2}={v}_{c}$ gives ${\rm{\Delta }}t\approx 1.2\,{t}_{\mathrm{ff}}({r}_{\max })$.

Assuming for the time being that the blob's cooling time remains nearly independent of gas pressure during its rise and fall, we can estimate the time of flight necessary for condensation by integrating $d\mathrm{ln}K/{dt}=-{t}_{\mathrm{cool}}^{-1}$ to obtain

$\begin{eqnarray}&&K\approx {K}_{1}{\left[1-\displaystyle \frac{6-3\lambda }{5}\displaystyle \frac{{\rm{\Delta }}t}{{t}_{\mathrm{cool}}({r}_{1})}\right]}^{5/(6-3\lambda )}.\end{eqnarray} \tag{ 7 }$

The condition for producing condensation is therefore

$\begin{eqnarray}&&{\rm{\Delta }}t\ \gtrsim \ \displaystyle \frac{5}{6-3\lambda }\,{t}_{\mathrm{cool}}({r}_{1}),\end{eqnarray} \tag{ 8 }$

or, equivalently,

$\begin{eqnarray}&&\displaystyle \frac{{v}_{1}}{{v}_{c}}\ \gtrsim \ \sqrt{2\ \mathrm{ln}\left[\displaystyle \frac{5}{6-3\lambda }\displaystyle \frac{{t}_{\mathrm{ff}}({r}_{\max })}{{\rm{\Delta }}t}{\left(\displaystyle \frac{{t}_{\mathrm{cool}}}{{t}_{\mathrm{ff}}}\right)}_{{r}_{1}}\right]}.\end{eqnarray} \tag{ 9 }$

The right-hand side of this second expression is generally close to unity, as long as ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ is not orders of magnitude greater than unity. Thus, the question of condensation depends primarily on the initial velocity of the blob's ballistic trajectory. Figure 2 shows the critical values of ${v}_{1}/{v}_{c}$ required for condensation, according to Equation (9), as functions of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ for $-1\leqslant \lambda \leqslant 0.5$. Notice that uplift velocities ${v}_{1}\approx {v}_{c}$ can induce condensation of blobs that initially have ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\approx 3$ but that ${v}_{1}\gtrsim 1.5{v}_{c}$ is needed for uplift to induce condensation of blobs that initially have ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\gtrsim 10$. This simple relationship is likely to be the biggest driver of condensation in more complex treatments of multiphase outflows (e.g., Thompson et al. 2016).

Zoom In Zoom Out Reset image size

Without gravity, the pressure P and entropy K of the background medium can be uniform. The contrast in entropy between the background and a perturbation with entropy $\tilde{K}$ is $\delta \mathrm{ln}K=\mathrm{ln}(\tilde{K}/K)$, and it evolves with time according to

$\begin{eqnarray}\displaystyle \frac{d}{{dt}}(\delta \mathrm{ln}K) & = & -\left(\displaystyle \frac{2}{3}\displaystyle \frac{\mu {m}_{p}}{k}\right)\delta \left(\displaystyle \frac{{ \mathcal L }}{T}\right)\\ & = & -\left(\displaystyle \frac{2}{3}\displaystyle \frac{\mu {m}_{p}}{k}\right){\left.\displaystyle \frac{\partial ({ \mathcal L }/T)}{\partial \mathrm{ln}K}\right|}_{A}\delta \mathrm{ln}K,\end{eqnarray} \tag{ 10 }$

where A is an arbitrary thermodynamic quantity complementary to K that remains constant during the process of condensation. Spatial gradients of ${ \mathcal L }$ do not appear in this equation because they must vanish in order for the background medium to remain uniform. The perturbation's amplitude therefore grows monotonically with time if

$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial ({ \mathcal L }/T)}{\partial \mathrm{ln}K}\right|}_{A}\ \lt \ 0,\end{eqnarray} \tag{ 11 }$

which is equivalent to the thermal instability condition originally derived by Balbus (1986). If the background medium is thermally balanced (${ \mathcal L }=0$), then this condition reduces to

$\begin{eqnarray}&&{\left.\displaystyle \frac{\partial { \mathcal L }}{\partial \mathrm{ln}K}\right|}_{A}\ \lt \ 0,\end{eqnarray} \tag{ 12 }$

which is equivalent to the thermal instability condition originally derived by Field (1965).

To obtain the timescale on which isobaric gas condenses, we set A = P, so that

$\begin{eqnarray}\displaystyle \frac{d}{{dt}}(\delta \mathrm{ln}K) & = & {\omega }_{\mathrm{ti}}(\delta \mathrm{ln}K),\\ {\omega }_{\mathrm{ti}} & \equiv & -\left(\displaystyle \frac{2}{3}\displaystyle \frac{\mu {m}_{p}}{k}\right){\left.\displaystyle \frac{\partial ({ \mathcal L }/T)}{\partial \mathrm{ln}K}\right|}_{P}.\end{eqnarray} \tag{ 13 }$

Then we separate ${ \mathcal L }$ into a cooling rate per unit mass C and a heating rate per unit mass H, so that ${ \mathcal L }=C-H$. If all of the cooling is radiative, then $C/T\propto {K}^{-(6-3\lambda )/5}{P}^{(1+2\lambda )/5}$ and

$\begin{eqnarray}&&{\omega }_{\mathrm{ti}}=\displaystyle \frac{6-3\lambda }{5\,{t}_{\mathrm{cool}}}\left[1+\left(\displaystyle \frac{5}{6-3\lambda }\right)\displaystyle \frac{H}{C}{\left.\displaystyle \frac{\partial \mathrm{ln}(H/T)}{\partial \mathrm{ln}K}\right|}_{P}\,\right].\end{eqnarray} \tag{ 14 }$

The second term inside the square brackets vanishes if H = 0. It also vanishes if heating per unit volume is constant, because in that case H/T is constant in an isobaric medium. Condensation therefore happens as long as $\lambda \lt 2$ and the factor in square brackets is positive, and it progresses exponentially on a timescale ${\omega }_{\mathrm{ti}}^{-1}$. Temporal variations in H have no effect on a perturbation's growth rate if both the perturbation and the ambient medium share the same value of H/T at each moment in time. However, small entropy perturbations in a medium with no heating (i.e., H = 0) do not get much of an opportunity to develop into nonlinear condensates, because the timescale for perturbation growth is then quite similar to the timescale for cooling of the background medium (e.g., Malagoli et al. 1987; Balbus & Soker 1989).

In the presence of gravity, the pressure of a static background medium cannot be uniform, and condensation couples interestingly with buoyancy (e.g., Defouw 1970; Cowie et al. 1980; Nulsen 1986; Malagoli et al. 1987; Balbus & Soker 1989; Loewenstein 1989). The strength of this coupling depends on the entropy gradient ${\rm{\nabla }}\mathrm{ln}K$ in the medium, because

$\begin{eqnarray}\displaystyle \frac{\partial }{\partial t}(\delta \mathrm{ln}K) & = & {\omega }_{\mathrm{ti}}(\delta \mathrm{ln}K)\\ & & +{\omega }_{\mathrm{sw}}(\delta \mathrm{ln}P)-\dot{{\boldsymbol{\xi }}}\cdot {\rm{\nabla }}\mathrm{ln}K,\end{eqnarray} \tag{ 15 }$

where ${\boldsymbol{\xi }}({\boldsymbol{r}},t)\equiv {\boldsymbol{r}}-{{\boldsymbol{r}}}_{0}({\boldsymbol{r}},t)$ represents the displacement of gas at location r and time t from its initial location ${{\boldsymbol{r}}}_{0}$ at time t = 0. In this expression, the quantity

$\begin{eqnarray}&&{\omega }_{\mathrm{sw}}\equiv -\left(\displaystyle \frac{2}{3}\displaystyle \frac{\mu {m}_{p}}{k}\right){\left.\displaystyle \frac{\partial ({ \mathcal L }/T)}{\partial \mathrm{ln}P}\right|}_{K,r}=-\displaystyle \frac{6+2\lambda }{5{t}_{\mathrm{cool}}}\end{eqnarray} \tag{ 16 }$

is the inverse of the timescale on which sound waves either grow or decay. The right-hand side of Equation (16) applies to a thermally balanced medium in which heating per unit volume is constant at each radius and implies that circumgalactic sound waves tend to be thermally stable, because ${\omega }_{\mathrm{sw}}\lt 0$ for essentially all environments of interest. Thermal instability consequently does not develop unless ${\omega }_{\mathrm{ti}}\gt 0$. However, subsequent motions resulting from changes in buoyancy can either increase or diminish the perturbation's (Eulerian) entropy contrast $\delta \mathrm{ln}K$, depending on the sign of $\dot{{\boldsymbol{\xi }}}\cdot {\rm{\nabla }}\mathrm{ln}K$.

We are primarily interested in systems in which ${ \mathcal L }$ depends explicitly on position, but it is illuminating to consider first the case in which ${ \mathcal L }$ depends on only the local thermodynamic state, as quantified by K and P. In that case, the thermal instability condition ${\omega }_{\mathrm{ti}}\gt 0$ cannot be satisfied in a hydrostatic medium that is convectively stable and in which sound waves are thermally stable, as emphasized by Balbus & Soker (1989). In such a medium, K increases with altitude while P decreases with altitude. Therefore, the condition ${ \mathcal L }(K,P)=0$ requires ${\omega }_{\mathrm{ti}}$ and ${\omega }_{\mathrm{sw}}$ to have the same sign. Otherwise, there cannot be a configuration that is simultaneously hydrostatic, convectively stable, and thermally balanced. Hence, a background configuration with all of three of these desirable stability characteristics does not exist unless ${\omega }_{\mathrm{ti}}\leqslant 0$.

An explicit dependence of ${ \mathcal L }$ on position ${\boldsymbol{r}}$ allows hydrostatic, convectively stable configurations to exist for a medium with ${\omega }_{\mathrm{ti}}\gt 0$ even if sound waves are thermally stable. In a spherically symmetric environment, a thermally balanced medium in which ${ \mathcal L }$ depends explicitly on position has

$\begin{eqnarray}&&{\alpha }_{K}{\left.\displaystyle \frac{\partial { \mathcal L }}{\partial \mathrm{ln}K}\right|}_{P,r}+\ {\alpha }_{P}{\left.\displaystyle \frac{\partial { \mathcal L }}{\partial \mathrm{ln}P}\right|}_{K,r}+\ {\left.\displaystyle \frac{\partial { \mathcal L }}{\partial \mathrm{ln}r}\right|}_{K,P}=0,\end{eqnarray} \tag{ 17 }$

where ${\alpha }_{K}\equiv {\boldsymbol{r}}\cdot {\rm{\nabla }}\mathrm{ln}K$ and ${\alpha }_{P}\equiv {\boldsymbol{r}}\cdot {\rm{\nabla }}\mathrm{ln}P$. The physical implications of this relationship are easiest to see if C depends only on K and P while H/T depends only on r, giving

$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}\,(H/T)}{d\mathrm{ln}r}=-({\alpha }_{K}\,{\omega }_{\mathrm{ti}}+{\alpha }_{P}\,{\omega }_{\mathrm{sw}})\,{t}_{\mathrm{cool}}.\end{eqnarray} \tag{ 18 }$

According to this equation, thermal instability can proceed in a thermally balanced medium satisfying the other three stability conditions as long as the decline in H/T with radius is sufficiently rapid.

From here on, we will focus on media that have ${\omega }_{\mathrm{ti}}\gt 0$ but in which sound waves are stable (${\omega }_{\mathrm{sw}}\lt 0$). Entropy perturbations can then be thermally unstable and undergo buoyancy-driven motions governed by the momentum equation

$\begin{eqnarray}(1+\delta \mathrm{ln}\rho )\ddot{{\boldsymbol{\xi }}} & = & \left[\displaystyle \frac{3}{5}(\delta \mathrm{ln}K)+\displaystyle \frac{2}{5}(\delta \mathrm{ln}P)\right]{\rm{\nabla }}\phi -\displaystyle \frac{P}{\rho }{\rm{\nabla }}(\delta \mathrm{ln}P).\end{eqnarray} \tag{ 19 }$

Here, ${\rm{\nabla }}\phi$ represents the gradient of the gravitational potential, and we have assumed ${\rm{\nabla }}P=-\rho {\rm{\nabla }}\phi$, so that the background state is in hydrostatic equilibrium. Differentiating Equation (15) with respect to time and applying Equation (19) after dropping the nonlinear $(\delta \mathrm{ln}\rho )\ddot{{\boldsymbol{\xi }}}$ term then gives

$\begin{eqnarray}&&\left[\,\displaystyle \frac{{\partial }^{2}}{\partial {t}^{2}}\ -\ {\omega }_{\mathrm{ti}}\displaystyle \frac{\partial }{\partial t}\ +\ {\omega }_{\mathrm{buoy}}^{2}\,\right](\delta \mathrm{ln}K)\\ &&\quad =\,\left[\,{\omega }_{\mathrm{sw}}\displaystyle \frac{\partial }{\partial t}-\,\displaystyle \frac{2}{3}{\omega }_{\mathrm{buoy}}^{2}+\,{\omega }_{\mathrm{buoy}}^{2}\displaystyle \frac{{c}_{s}^{2}}{g}\,\displaystyle \frac{\partial }{\partial r}\right](\delta \mathrm{ln}P),\end{eqnarray} \tag{ 20 }$

where ${\omega }_{\mathrm{buoy}}\equiv {[{(3/5){\rm{\nabla }}\phi \cdot {\rm{\nabla }}\mathrm{ln}K]}^{1/2}=(6{\alpha }_{K}/5)}^{1/2}{t}_{\mathrm{ff}}^{-1}$ is the usual Brunt–Väisälä frequency for buoyant oscillations and ${c}_{s}={(5P/3\rho )}^{1/2}$ is the adiabatic sound speed. Note that this equation has been obtained without using either the continuity equation or the Boussinesq approximation and relates internal gravity waves (left-hand side) to the corresponding pressure disturbances (right-hand side).

If the right-hand side of the equation is negligible, then an entropy perturbation can move through the background medium without dissipative resistance, and the qualitative behavior of the perturbation's (Eulerian) entropy contrast $\delta \mathrm{ln}K$ is obvious. Entropy perturbations in a medium with ${\omega }_{\mathrm{ti}}\gg {\omega }_{\mathrm{buoy}}\gt 0$ grow monotonically with time, while those in a medium with ${\omega }_{\mathrm{buoy}}\gg {\omega }_{\mathrm{ti}}\gt 0$ oscillate with continually increasing amplitude. The oscillations grow because of a phase difference between ${\xi }_{r}\equiv ({\boldsymbol{r}}\cdot {\boldsymbol{\xi }})/r$, which is the perturbation's displacement in the radial direction and the radial restoring force, which is proportional to $\delta \mathrm{ln}K$. During most of each oscillation cycle, the product of ${\xi }_{r}$ and $\delta \mathrm{ln}K$ is negative, meaning that buoyancy forces are usually pushing the perturbation back toward the midpoint of the oscillations. However, $\delta \mathrm{ln}K\gt 0$ at the beginning of each upward excursion above the midpoint, and $\delta \mathrm{ln}K\lt 0$ at the beginning of each downward excursion below the midpoint. The resulting buoyancy impulses cause the amplitudes of both ${\boldsymbol{\xi }}$ and $\delta \mathrm{ln}K$ to increase with time, as originally found by Defouw (1970).

Substituting $\delta \mathrm{ln}K\propto \,{e}^{-i\omega t}$ into Equation (20) with the right-hand side set to zero gives the frequency solutions

$\begin{eqnarray}&&\omega \ =\ {\omega }_{\pm }\ \equiv \ \displaystyle \frac{i{\omega }_{\mathrm{ti}}}{2}\pm {\left({\omega }_{\mathrm{buoy}}^{2}-\displaystyle \frac{{\omega }_{\mathrm{ti}}^{2}}{4}\right)}^{1/2}.\end{eqnarray} \tag{ 21 }$

Low-entropy perturbations in a such a medium progress steadily toward nonlinear condensation if ${\omega }_{\mathrm{buoy}}\lt {\omega }_{\mathrm{ti}}/2$, which corresponds to

$\begin{eqnarray}&&{\alpha }_{K}\lt \displaystyle \frac{3{(2-\lambda )}^{2}}{40}{\left(\displaystyle \frac{{t}_{\mathrm{ff}}}{{t}_{\mathrm{cool}}}\right)}^{2},\end{eqnarray} \tag{ 22 }$

when expressed in terms of ${t}_{\mathrm{cool}}$, ${t}_{\mathrm{ff}}$, and ${\alpha }_{K}$. However, the qualitative behavior of the oscillatory entropy perturbations that arise when ${\omega }_{\mathrm{buoy}}\gg {\omega }_{\mathrm{ti}}$ depends critically on how one treats the $\delta \mathrm{ln}P$ side of the equation. We will now proceed with a careful look at those $\delta \mathrm{ln}P$ terms, in order to clarify some seemingly contradictory claims in the astrophysical literature on local thermal instability.

3.3.1. The Bobbing-blob Approximation

One intuitively appealing approach is to idealize an entropy perturbation as a discrete, isolated blob in which all of the gas shares the same displacement vector (e.g., Cowie et al. 1980; Nulsen 1986; Loewenstein 1989). In a medium with ${\omega }_{\mathrm{buoy}}\gg {\omega }_{\mathrm{ti}}$, the vertical bobbing of the blob necessarily excites internal gravity waves that propagate away from the blob. This transfer of kinetic energy from the blob to its environment implies that the blob's motion must experience some resistance in addition to the background pressure gradient responsible for maintaining the oscillations. Adding the necessary resistance by replacing the pressure and density perturbation terms in the equation of motion with a heuristic dissipation term, so that

$\begin{eqnarray}&&\ddot{{\boldsymbol{\xi }}}=\displaystyle \frac{3}{5}(\delta \mathrm{ln}K){\rm{\nabla }}\phi \,-\,{\omega }_{{\rm{D}}}\dot{{\boldsymbol{\xi }}},\end{eqnarray} \tag{ 23 }$

qualitatively accounts for those kinetic energy losses by dissipating the blob's kinetic energy on the timescale ${\omega }_{{\rm{D}}}^{-1}$. This replacement results in the frequency solutions

$\begin{eqnarray}&&\omega ={\omega }_{\pm }\equiv \displaystyle \frac{i({\omega }_{\mathrm{ti}}-{\omega }_{{\rm{D}}})}{2}\pm \ {\left[{\omega }_{\mathrm{buoy}}^{2}-\displaystyle \frac{{({\omega }_{\mathrm{ti}}+{\omega }_{{\rm{D}}})}^{2}}{4}\right]}^{1/2},\end{eqnarray} \tag{ 24 }$

which show that the blob's oscillations decay instead of growing if the timescale for dissipating kinetic energy is shorter than the thermal timescale for pumping the oscillations.

Associating the dissipation term with hydrodynamical drag gives ${\omega }_{{\rm{D}}}\sim | \dot{{\boldsymbol{\xi }}}| /{\rm{\Delta }}r$, where ${\rm{\Delta }}r$ is the blob's thickness in the vertical direction. Given this scaling for ${\omega }_{{\rm{D}}}$, thermal pumping will increase the amplitudes of small oscillations until they saturate at $\delta \mathrm{ln}K\sim {\alpha }_{K}({\omega }_{\mathrm{ti}}/{\omega }_{\mathrm{buoy}})({\rm{\Delta }}r/r)$. The largest (Eulerian) density fluctuations, $\delta \mathrm{ln}\rho \sim {\alpha }_{K}^{1/2}({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})$, are therefore obtained for ${\rm{\Delta }}r\sim r$. This analytical estimate of the saturation threshold agrees well with numerical simulations of oscillatory thermal instability in thermally balanced media (e.g., McCourt et al. 2012; Meece et al. 2015). However, the route through which it was obtained indicates that the saturation process is inherently nonlinear. Furthermore, the bobbing-blob approximation ignores what happens to the propagating internal gravity waves stimulated by the blob's oscillations.

3.3.2. Plane-wave Perturbations

Alternatively, one can treat the perturbations as Eulerian plane waves (e.g., Malagoli et al. 1987; Balbus 1988; Balbus & Soker 1989), but then some care needs to be taken with the physical interpretation. Making the usual plane-wave assumption that $\delta \mathrm{ln}K$, $\delta \mathrm{ln}P$, and ${\boldsymbol{\xi }}$ are all $\propto \exp [i({\boldsymbol{k}}\cdot {\boldsymbol{r}}-\omega t)]$ and defining ${ \mathcal R }\equiv \delta \mathrm{ln}P/\delta \mathrm{ln}K$ allows us to convert Equation (20) to the dispersion relation

$\begin{eqnarray}{\omega }^{2} & - & i\,({\omega }_{\mathrm{ti}}+\,{ \mathcal R }{\omega }_{\mathrm{sw}})\,\omega \\ & - & \left(1+\displaystyle \frac{2}{3}{ \mathcal R }-i{ \mathcal R }\displaystyle \frac{{c}_{s}^{2}{k}_{r}}{g}\right){\omega }_{\mathrm{buoy}}^{2}=0.\end{eqnarray} \tag{ 25 }$

In order to account for the influence of pressure perturbations on this dispersion relation through the ratio ${ \mathcal R }$, we need to bring in the continuity equation for the perturbation, which can be expressed as

$\begin{eqnarray}&&\displaystyle \frac{5}{3}({\rm{\nabla }}\cdot {\boldsymbol{\xi }})\ =\ {\rm{\Delta }}\mathrm{ln}K-{\rm{\Delta }}\mathrm{ln}P\end{eqnarray} \tag{ 26 }$

when written in terms of the Lagrangian perturbation amplitudes ${\rm{\Delta }}\mathrm{ln}P=\delta \mathrm{ln}P+{\alpha }_{P}({\xi }_{r}/r)$ and ${\rm{\Delta }}\mathrm{ln}K\equiv \delta \mathrm{ln}K+{\alpha }_{K}({\xi }_{r}/r)$. Using the plane-wave assumption to simplify the tangential momentum equation, one can write the divergence of the displacement field as

$\begin{eqnarray}&&{\rm{\nabla }}\cdot {\boldsymbol{\xi }}=({{ik}}_{r}r+{\rm{\nabla }}\cdot {\boldsymbol{r}})\displaystyle \frac{{\xi }_{r}}{r}-\displaystyle \frac{3{c}_{s}^{2}{k}_{\perp }^{2}}{5{\omega }^{2}}\delta \mathrm{ln}P,\end{eqnarray} \tag{ 27 }$

where k_r and ${k}_{\perp }$ are the radial and tangential wavenumbers, respectively. The coordinate divergence term is ${\rm{\nabla }}\cdot {\boldsymbol{r}}=2$ in a spherical environment but can be ignored in a plane-parallel approximation. Combining Equations (26) and (27) and recognizing that ${\alpha }_{P}=-5{gr}/3{c}_{s}^{2}$ then leads to

$\begin{eqnarray}{ \mathcal R } & = & \displaystyle \frac{{\omega }^{2}}{{\omega }^{2}-{c}_{s}^{2}{k}_{\perp }^{2}}\left[\displaystyle \frac{5}{3{\alpha }_{K}}\left({{ik}}_{r}r+{\rm{\nabla }}\cdot {\boldsymbol{r}}-\displaystyle \frac{{gr}}{{c}_{s}^{2}}\right)\right.\\ & & \times \left.\left(1-\displaystyle \frac{{\rm{\Delta }}\mathrm{ln}K}{\delta \mathrm{ln}K}\right)+\displaystyle \frac{{\rm{\Delta }}\mathrm{ln}K}{\delta \mathrm{ln}K}\right].\end{eqnarray} \tag{ 28 }$

Applying the plane-wave assumption to Equation (15) gives the ratio of Lagrangian to Eulerian entropy perturbation amplitudes:

$\begin{eqnarray}&&\displaystyle \frac{{\rm{\Delta }}\mathrm{ln}K}{\delta \mathrm{ln}K}=i\left(\displaystyle \frac{{\omega }_{\mathrm{ti}}+{ \mathcal R }{\omega }_{\mathrm{sw}}}{\omega }\right).\end{eqnarray} \tag{ 29 }$

The dispersion relation for perturbations with short radial wavelengths ($| {k}_{r}r| \gg 1$) therefore reduces to

$\begin{eqnarray}{\omega }^{2} & - & i({\omega }_{\mathrm{ti}}+{ \mathcal R }{\omega }_{\mathrm{sw}})\omega \\ & - & \left(1+\displaystyle \frac{2}{3}{ \mathcal R }\right)\left(\displaystyle \frac{{c}_{s}^{2}{k}_{\perp }^{2}-{\omega }^{2}}{{c}_{s}^{2}{k}^{2}-{\omega }^{2}}\right){\omega }_{\mathrm{buoy}}^{2}=0,\end{eqnarray} \tag{ 30 }$

where ${k}^{2}\equiv {k}_{r}^{2}+{k}_{\perp }^{2}$.

The literature on thermal instability in circumgalactic media has generally focused on short-wavelength modes with ${\omega }^{2}\ll {c}_{s}^{2}{k}_{\perp }^{2}/| {k}_{r}r|$, which are nearly pure internal gravity waves. In that case, the amplitude ratio ${ \mathcal R }=\delta \mathrm{ln}P/\delta \mathrm{ln}K$ is small, and the dispersion relation reduces further to

$\begin{eqnarray}&&{\omega }^{2}\,-\,i\,{\omega }_{\mathrm{ti}}\omega \,-\displaystyle \frac{{k}_{\perp }^{2}}{{k}^{2}}{\omega }_{\mathrm{buoy}}^{2}=0\end{eqnarray} \tag{ 31 }$

(e.g., Malagoli et al. 1987; Binney et al. 2009). For these modes, Eulerian pressure perturbations influence the gravity-wave oscillator only through the pressure-gradient term in Equation (20). The resulting ${k}_{\perp }^{2}/{k}^{2}$ factor slows the oscillations because the gravitational force correcting a given radial displacement meets with greater inertial resistance from transverse displacements as ${k}_{\perp }^{2}/{k}^{2}$ declines (e.g., Binney et al. 2009). As a consequence, the critical value of ${\alpha }_{K}$ separating monotonic perturbation growth from oscillatory growth for a given ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ ratio becomes larger than one would obtain by setting $\delta \mathrm{ln}P=0$.

Some treatments in the literature (e.g., Balbus & Soker 1989; McCourt et al. 2012) arrive at a similar dispersion relation by setting ${\rm{\nabla }}\cdot {\boldsymbol{\xi }}=0$. With that restriction, we find

$\begin{eqnarray}&&{ \mathcal R }=\displaystyle \frac{{\alpha }_{P}}{{\alpha }_{K}}\left(1-\displaystyle \frac{{\rm{\Delta }}\mathrm{ln}K}{\delta \mathrm{ln}K}\right)+\displaystyle \frac{{\rm{\Delta }}\mathrm{ln}K}{\delta \mathrm{ln}K},\end{eqnarray} \tag{ 32 }$

and obtain the dispersion relation

$\begin{eqnarray}&&{\omega }^{2}-i({\omega }_{\mathrm{ti}}+\,{ \mathcal R }{\omega }_{\mathrm{sw}})\,\omega -\left(1+\displaystyle \frac{2}{3}{ \mathcal R }\right)\\ &&\quad \,\times \,\left[\displaystyle \frac{{k}_{\perp }^{2}{r}^{2}}{{k}^{2}{r}^{2}+{{ik}}_{r}r({\rm{\nabla }}\cdot {\boldsymbol{r}})}\right]{\omega }_{\mathrm{buoy}}^{2}=0\end{eqnarray} \tag{ 33 }$

without having to assume $| {k}_{r}r| \gg 1$. If the ${{ik}}_{r}r({\rm{\nabla }}\cdot {\boldsymbol{r}})$ term can be ignored, then this relation would reduce to Equation (31) if ${ \mathcal R }$ were small. However, $| { \mathcal R }|$ is typically of order unity for divergence-free modes, indicating that they are qualitatively different from the internal gravity waves discussed in the previous paragraph. In fact, the restriction ${\rm{\nabla }}\cdot {\boldsymbol{\xi }}=0$ implies that the radial wavenumber for such modes is imaginary, with

$\begin{eqnarray}&&{k}_{r}r=i\left[{k}_{\perp }^{2}{r}^{2}\left(\displaystyle \frac{{c}_{s}^{2}}{{gr}}\right)\displaystyle \frac{{\omega }_{\mathrm{buoy}}^{2}}{{\omega }^{2}}+({\rm{\nabla }}\cdot {\boldsymbol{r}})\right]\end{eqnarray} \tag{ 34 }$

in the adiabatic limit. These divergence-free modes are therefore evanescent in the radial direction. Furthermore, the associated growth rate of thermal instability is ${\omega }_{\mathrm{ti}}+({\alpha }_{P}/{\alpha }_{K}){\omega }_{\mathrm{sw}}$, which exceeds ${\omega }_{\mathrm{ti}}$ in media with stable sound waves (i.e., ${\omega }_{\mathrm{sw}}\lt 0$), because ${\alpha }_{P}/{\alpha }_{K}\lt 0$.

3.3.3. Onset of Nonlinearity

Prior analyses of circumgalactic thermal instability did not focus on the consequences of modes in which ${\omega }_{\mathrm{sw}}$ is important, but here we would like to understand why the growth of oscillatory thermal instability saturates with $\delta \mathrm{ln}\rho \,\sim {\alpha }_{K}^{1/2}({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})$, as predicted by the heuristic bobbing-blob model and as observed in numerical simulations (McCourt et al. 2012; Meece et al. 2015). If only low-${ \mathcal R }$ or divergence-free modes are permitted, then there is no channel (other than viscosity) for the dissipation of the kinetic energy that oscillatory thermal instability introduces into the system. Nonlinear mode coupling can exchange energy among those gravity-wave modes but does not by itself limit the overall growth of a set of modes that are all thermally unstable. In order to achieve a steady state in which $\delta \mathrm{ln}\rho$ depends on ${t}_{\mathrm{cool}}$, nonlinear mode coupling must instead channel the accumulating gravity-wave energy into acoustic modes, which damp through radiative losses on a timescale $\sim | {\omega }_{\mathrm{sw}}^{-1}| \sim {t}_{\mathrm{cool}}$.

In lieu of presenting a complete solution to that mode-coupling problem in this paper, we would instead like to briefly illustrate how the gravity-wave saturation scale $\delta \mathrm{ln}\rho \,\sim {\alpha }_{K}^{1/2}({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})$ arises from nonlinearity. We are interested in knowing how mode coupling changes the amplitudes A_i of the radial-displacement eigenmodes ${\xi }_{i}={A}_{i}\exp [i({{\boldsymbol{k}}}_{i}\cdot {\boldsymbol{r}}-{\omega }_{i}t)]$, for which ${\omega }_{i}({{\boldsymbol{k}}}_{i})$ satisfies Equations (25), (28), and (29). The full radial momentum equation for a superposition of those modes can be expressed as

$\begin{eqnarray}&&(1+\delta \mathrm{ln}\rho )\left(\displaystyle \sum _{i}{\ddot{\xi }}_{i}\right)\ =\ -\displaystyle \sum _{j}{\omega }_{j}^{2}{\xi }_{j},\end{eqnarray} \tag{ 35 }$

and the amplitudes A_i remain constant if the nonlinear terms are ignored. After canceling the linear terms, we can rewrite this equation in a form that relates the rates of amplitude change to those nonlinear terms,

$\begin{eqnarray}&&\displaystyle \sum _{i}\left[\left(\displaystyle \frac{{\ddot{A}}_{i}-i2{\omega }_{i}{\dot{A}}_{i}}{{A}_{i}}\right){\xi }_{i}\right]\\ &&\quad \quad =\,-\,\left(\displaystyle \sum _{j}{\omega }_{j}^{2}{\xi }_{j}\right)\left[\displaystyle \sum _{n=1}^{\infty }{(-\delta \mathrm{ln}\rho )}^{n}\right]\\ &&\quad \quad \approx \,-\left(\displaystyle \sum _{j}{\omega }_{j}^{2}{\xi }_{j}\right)\left[\displaystyle \sum _{m}{(\delta \mathrm{ln}\rho )}_{m}\right],\end{eqnarray} \tag{ 36 }$

where ${(\delta \mathrm{ln}\rho )}_{m}$ represents the oscillating density perturbation corresponding to eigenmode ${\xi }_{m}$, and we have kept only the lowest-order nonlinear terms on the right-hand side. Resonant energy transfer can therefore happen among triads of modes that satisfy ${{\boldsymbol{k}}}_{i}={{\boldsymbol{k}}}_{j}+{{\boldsymbol{k}}}_{m}$ and ${\omega }_{i}={\omega }_{j}+{\omega }_{m}$.

Here, we are most interested in coupling between a typical internal gravity wave of wavenumber k and frequency ${\omega }_{g}\approx \,({k}_{\perp }/k){\omega }_{\mathrm{buoy}}$, and a pair of sound waves with frequencies $\sim {c}_{s}k$ that differ by ${\omega }_{g}$ so that their beat pattern resonates with the gravity wave. In a purely adiabatic system, this triad of waves would eventually come into energy equipartition with ${A}_{g}{\omega }_{g}\sim {A}_{s}{c}_{s}k$, where A_g and A_s are the gravity-wave and sound-wave displacement amplitudes, respectively. On the right-hand side of Equation (36), terms corresponding to pairs of sound waves of appropriate wavelength have magnitudes $\sim {({A}_{s}{c}_{s}k)}^{2}k$ because ${(\delta \mathrm{ln}\rho )}_{m}\sim {A}_{s}{k}_{m}$. In a system in equipartition, these terms are $\sim {({A}_{g}{\omega }_{g})}^{2}k$, implying that mode coupling transfers energy from sound waves into gravity waves and back again at a rate $\sim {A}_{g}{\omega }_{g}k$.

The system we are considering is not adiabatic. Instead, gravity waves grow in amplitude on a timescale $\sim {\omega }_{\mathrm{ti}}^{-1}$ until their amplitudes are large enough for mode coupling to channel their energy into dissipative sound waves on the same timescale. This happens when

$\begin{eqnarray}\displaystyle \frac{{A}_{g}}{r} & \sim & \displaystyle \frac{{\omega }_{\mathrm{ti}}}{{\omega }_{\mathrm{buoy}}}\displaystyle \frac{1}{{kr}},\\ {(\delta \mathrm{ln}\rho )}_{g} & \sim & \delta \mathrm{ln}K\ \sim \ {\alpha }_{K}\displaystyle \frac{{A}_{g}}{r}\\ & \sim & {\alpha }_{K}\displaystyle \frac{{\omega }_{\mathrm{ti}}}{{\omega }_{\mathrm{buoy}}}\displaystyle \frac{1}{{kr}}\ \sim \ {\alpha }_{K}^{1/2}\displaystyle \frac{{t}_{\mathrm{ff}}}{{t}_{\mathrm{cool}}}\displaystyle \frac{1}{{kr}},\end{eqnarray} \tag{ 37 }$

where ${(\delta \mathrm{ln}\rho )}_{g}$ is the density perturbation scale associated with gravity-wave oscillations.

The saturation amplitude derived from nonlinear mode coupling considerations in Section 3.3.3 closely resembles the one derived from the heuristic bobbing-blob model in Section 3.3.1. In the bobbing-blob case, the timescale for the energy transfer from the blob to the surrounding medium is $\sim {({A}_{g}{\omega }_{\mathrm{buoy}})}^{-1}{\rm{\Delta }}r$, and energy is permanently lost if it propagates away from the blob without contributing to the restoring force on it. Analogously, the timescale for energy transfer from internal gravity waves to damped acoustic modes is $\sim {({A}_{g}{\omega }_{\mathrm{buoy}})}^{-1}\,{k}^{-1}$, regardless of the mechanism that dissipates sound-wave energy. In both cases, the saturation amplitudes are largest for inhomogeneities with ${kr}\sim 1$.

Saturation of thermally unstable gravity-wave oscillations therefore has important consequences for condensation in the circumgalactic medium. In an otherwise static medium with ${\alpha }_{K}\gg {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$, energy transfer from buoyant oscillations into dissipative modes prevents a gas blob with a slightly lower entropy than its surroundings from condensing. At first, the blob descends faster than it can radiate away its thermal energy and accelerates until it passes through a layer of equivalent entropy. After it overshoots that layer, it begins to bob, and its subsequent behavior depends on the amplitude $| \delta \mathrm{ln}K|$ of its Eulerian entropy contrast. If that amplitude is less than the saturation amplitude, then thermal pumping will cause it to grow to the saturation amplitude. But if $| \delta \mathrm{ln}K|$ is greater than the saturation amplitude, then coupling of buoyant oscillations to acoustic modes causes $| \delta \mathrm{ln}K|$ to decay until it reaches the saturation amplitude. Hereafter, we will refer to such a decline in perturbation amplitude as buoyancy damping, because it results from the damping of buoyant oscillations in a gravitationally stratified medium.

In contrast, buoyancy does not inhibit condensation in a medium with ${\alpha }_{K}\ll {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$, because there are no oscillations. Instead, low-entropy blobs proceed steadily toward becoming nonlinear condensates. In fact, the onset of nonlinearity can actually assist perturbation growth, because it inhibits the descent of low-entropy gas toward layers of equivalent entropy (e.g., Nulsen 1986). Therefore, the condition

$\begin{eqnarray}&&{\alpha }_{K}\sim {\left(\displaystyle \frac{{t}_{\mathrm{ff}}}{{t}_{\mathrm{cool}}}\right)}^{2}\end{eqnarray} \tag{ 38 }$

corresponds to a critical entropy-profile slope separating media in which linear thermal instability leads to condensation from media in which linear thermal instability saturates at relatively small amplitudes. ⁵

At the margin between these regimes, where ${\omega }_{\mathrm{buoy}}\sim {\omega }_{\mathrm{ti}}$, thermal instability itself starts to reduce the frequency of gravity-wave oscillations, which slows the rate of energy transfer into sound waves. In a system with ${\alpha }_{K}\gtrsim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$, the saturation scale for thermally unstable gravity-wave oscillations with ${kr}\sim 1$ becomes

$\begin{eqnarray}\displaystyle \frac{{A}_{g}}{r} & \sim & \displaystyle \frac{{\omega }_{\mathrm{ti}}}{{\omega }_{\mathrm{buoy}}}{\left(1-\displaystyle \frac{{\omega }_{\mathrm{ti}}^{2}}{4{\omega }_{\mathrm{buoy}}^{2}}\right)}^{-1/2},\\ {(\delta \mathrm{ln}\rho )}_{g} & \sim & {\alpha }_{K}\displaystyle \frac{{A}_{g}}{r}\ \gtrsim \ \displaystyle \frac{{t}_{\mathrm{ff}}^{2}}{{t}_{\mathrm{cool}}^{2}}{\left(1-\displaystyle \frac{{\omega }_{\mathrm{ti}}^{2}}{4{\omega }_{\mathrm{buoy}}^{2}}\right)}^{-1/2}.\end{eqnarray} \tag{ 39 }$

In other words, there is a lower limit of ${(\delta \mathrm{ln}\rho )}_{g}\sim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$ on the saturation amplitude of gravity-wave perturbations, and the limiting amplitude approaches unity as ${\omega }_{\mathrm{buoy}}$ approaches ${\omega }_{\mathrm{ti}}/2$ from above. A reduction of the entropy gradient to ${\alpha }_{K}\sim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$ consequently eliminates buoyancy damping and allows thermal instability to proceed to condensation. ⁶

McCourt et al. (2012) drew attention to the dependence of saturation amplitude on ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ using simulations of a plane-parallel medium in which heating balances cooling within each equipotential layer. They showed that small entropy perturbations introduced into an initially isothermal medium grow to an Eulerian entropy contrast $\sim {\alpha }_{K}({\omega }_{\mathrm{ti}}/{\omega }_{\mathrm{buoy}})$ before saturating. The interpretation they offered for this saturation was that nonlinear mode coupling leads to dissipation when ${\xi }_{r}/r\sim {\omega }_{\mathrm{ti}}/{\omega }_{\mathrm{buoy}}$, but they were not specific about the nature of the dissipation or the mechanism through which mode coupling disrupts the growth of gravity-wave oscillations. They concluded that the criterion ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\lesssim 1$ was a necessary condition for condensation in a plane-parallel medium.

Meece et al. (2015) performed similar simulations showing that the condensation criterion is more complex, in that condensation is always possible near the midplane, even in systems with ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\gg 1$, if enough time is allowed (see also Choudhury & Sharma 2016). Figure 3 shows some previously unpublished results from the Meece et al. effort. Perturbation growth far from the midplane saturates as in McCourt et al. (2012) but is not fully suppressed near the midplane. Instead, it is significantly delayed for large values of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\gg 1$ and eventually generates crosstalk that affects perturbation growth in higher-lying layers, for reasons we will now proceed to analyze.

Zoom In Zoom Out Reset image size

Condensation is always possible near the midplane of a macrophysically balanced system in a plane-parallel potential well because buoyancy damping is ineffective there as long as the gas density does not become infinite at the midplane. In the bobbing-blob approximation of Section 3.3.1, the boundary between monotonic growth and oscillatory thermal instability is at ${\omega }_{D}={\omega }_{\mathrm{ti}}$, implying that the condition for condensation is

$\begin{eqnarray}&&{\alpha }_{K}\lt \displaystyle \frac{3{(2-\lambda )}^{2}}{10}{\left(\displaystyle \frac{{t}_{\mathrm{ff}}}{{t}_{\mathrm{cool}}}\right)}^{2}.\end{eqnarray} \tag{ 40 }$

In order to determine the zone around the midplane in which condensation can proceed, one can adopt a linear approximation to the potential, in which the freefall time ${t}_{\mathrm{ff},0}$ is constant with height. The pressure and entropy at the midplane are P₀ and K₀, respectively, giving a midplane density ${\rho }_{0}\,=\mu {m}_{p}{(n/{n}_{e})}^{2/5}{({P}_{0}/{K}_{0})}^{3/5}$. To lowest order in height z above the midplane, the pressure profile is $P{(z)={P}_{0}[1-(z/{z}_{P})}^{2}]$, where ${z}_{P}\equiv {t}_{\mathrm{ff},0}{({P}_{0}/{\rho }_{0})}^{1/2}$, and the linearized entropy profile can be expressed as $K(z)={K}_{0}[1+(z/{z}_{K})]$. The condensation condition in this zone is therefore

$\begin{eqnarray}\displaystyle \frac{z}{{z}_{K}} & \lesssim & \displaystyle \frac{3{(2-\lambda )}^{2}}{10}{\left[\displaystyle \frac{{t}_{\mathrm{ff},0}}{{t}_{\mathrm{cool},0}}\right]}^{2}{\left[1+\left(\displaystyle \frac{z}{{z}_{K}}\right)\right]}^{6(2-\lambda )/5}\\ & & \times {\left[1-{\left(\displaystyle \frac{z}{{z}_{P}}\right)}^{2}\right]}^{2(1+2\lambda )/5},\end{eqnarray} \tag{ 41 }$

where ${t}_{\mathrm{cool},0}\equiv {t}_{\mathrm{cool}}({K}_{0},{P}_{0})$. In other words, the thickness of the midplane region that is unstable to condensation is approximately ${({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$ times the entropy scale height in media with ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\gg 1$.

At first glance, the susceptibility to condensation of a system in which ${\alpha }_{K}\lesssim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$ would appear to be in conflict with the plane-parallel simulations of McCourt et al. (2012) that start with isentropic initial conditions (${\alpha }_{K}=0$) and fail to produce nonlinear density contrasts. In those simulations, the fractional fluctuations in entropy and density at one pressure scale height above the midplane saturate at $\sim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$ and remain at this amplitude for many cooling times. The reason for saturation of perturbation growth in what would seem to be an isentropic medium is that buoyant migration of the initial perturbations causes an entropy gradient to develop. As lower-entropy fluctuations sink and higher-entropy fluctuations rise, the resulting radial sorting of entropy perturbations increases ${\alpha }_{K}$. Eventually, the entropy gradient exceeds ${\alpha }_{K}\sim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$, and the condensation condition is no longer satisfied. The system then transitions to the bobbing-blobs regime, in which (Eulerian) entropy fluctuations saturate at a fractional amplitude $\sim {\alpha }_{K}({\omega }_{\mathrm{ti}}/{\omega }_{\mathrm{buoy}})$, which corresponds to $\sim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$ for ${\alpha }_{K}\sim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$, as shown in Section 3.4.

McCourt et al. (2012) recognized that the development of an entropy gradient through convective sorting of the initial perturbations was responsible for suppressing perturbation growth away from the midplane but did not fully explore the consequences for stability close to the midplane. Instead, they suppressed both heating and cooling near the midplane. In contrast, Meece et al. (2015) chose to allow heating and cooling to proceed near the midplane, and those simulations ended up displaying qualitative behavior distinctly different from the initially isentropic McCourt et al. (2012) simulations, even though they started with nearly isentropic initial conditions within a pressure scale height of the midplane.

Midplane condensation always occurs in the Meece et al. (2015) simulations, no matter what the initial value of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$, because there is always a zone around the midplane in which ${\alpha }_{K}\lesssim {({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})}^{2}$. However, buoyancy damping delays condensation by a time interval depending on ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$. Without gravity and buoyancy damping, the initial entropy perturbations of Meece et al. (2015) would experience monotonic exponential growth. Nearly exponential growth is indeed seen in simulations with ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\lesssim 3$, but larger values of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ lead to the saturation of perturbation growth at progressively smaller amplitudes (see Figure 3). Buoyancy damping suppresses perturbation growth away from the midplane but succeeds only in delaying it near the midplane. Delays occur because buoyancy damping limits the initial amplitudes of the seed perturbations that enter the isentropic zone and ultimately condense. As a result, the time required for condensation at the midplane is a larger multiple of ${t}_{\mathrm{cool}}$ in systems with larger values of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$.

Once condensation begins at the midplane in the Meece et al. (2015) simulations, the idealized algorithm for enforcing thermal balance within each layer ends up stimulating perturbation growth at greater altitudes. As condensates form near the midplane, the requirement for heating to balance cooling in each layer causes compensating hot bubbles to form. Those hot bubbles then rise out of the condensation zone. When they reach higher layers, the bubbles catalyze additional condensation, because regions that were formerly at the mean entropy for their layer, and therefore close to thermal balance, suddenly find themselves significantly below the new mean entropy level. The net cooling rate in lower-entropy regions of the layer must therefore increase in order to satisfy the thermal balance constraint. As a result, convective crosstalk between layers causes nonlinear perturbation growth to propagate beyond the thermally unstable midplane zone. The same effect was also found in simulations by Choudhury & Sharma (2016).

Does this crosstalk phenomenon, which arises in Meece et al. (2015) and Choudhury & Sharma (2016) from a highly idealized representation of global thermal balance, have a qualitative counterpart in more realistic treatments of feedback? We suggest that the answer to this question may be yes, because the outflow of feedback energy from a central source in a real galactic system is likely to propagate primarily through the lowest-density, highest-entropy ambient gas as it follows the path of least resistance to greater altitudes. Unless some mechanism (such as thermal conduction) exists for preferentially channeling feedback energy into denser, cooler gas clumps, then lower-entropy regions will inevitably cool and condense while feedback energy flows outward around them, even if heating and cooling are in global balance.

This phenomenon illustrates how centrally injected thermal feedback renders the surrounding region vulnerable to condensation: convection eliminates the central entropy gradient and therefore shuts off buoyancy damping of perturbation growth. Consequently, centrally injected thermal feedback actually promotes precipitation rather than suppressing it. This is why centrally injected thermal feedback fails so miserably to prevent condensation in the cluster-core simulation of Meece et al. (2017) that relies on pure thermal feedback, as well as many other such simulations that rely on central thermal feedback alone. With that thought in mind, we now turn our attention from local thermal instability to its implications for global thermal balance.

From a thermodynamic perspective, development of a multiphase medium under conditions of global thermal balance is essentially the opposite of thermal mixing in a system with constant total energy. For example, consider a system in which a small fraction of the gas condenses by radiating an amount of heat energy $\delta Q$. To maintain strict thermal balance, the rest of the gas in the system must gain the same amount of heat energy through feedback. If this condensation and feedback process proceeds at constant pressure, then the volume of the overall system does not change, but its overall entropy must decrease, because hotter gas gains heat energy while colder gas radiates it away. This inhomogeneous condensation process is prohibited in a closed system by the second law of thermodynamics. However, the supply of free energy that is added to balance cooling in a globally balanced system enables it to "unmix" via condensation. In this context, feedback is not only necessary to keep the system in overall thermal balance, it is also required to transform the original single-phase medium into a multiphase medium via thermal instability.

Gravity can suppress "unmixing" by causing buoyancy damping. The simulation results discussed in Section 4 show that buoyancy damping in a gravitationally stratified medium causes thermal instability to saturate in layers with ${\alpha }_{K}\gg {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$ and can substantially delay condensation even if ${\alpha }_{K}\ll {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$. Extrapolating these findings to a cluster-core environment suggests that a thermally balanced cluster core with ${t}_{\mathrm{cool}}\lesssim 1\,\mathrm{Gyr}$ and $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\sim 10$ may be able to resist condensation for $\gtrsim 5\,\mathrm{Gyr}$.

The idealized cluster-core simulations of Sharma et al. (2012b), which strictly maintain thermal balance in each equipotential layer, appear to corroborate this suggestion. In their fiducial simulation (run C10), the core configuration begins with central entropy $\approx 10\,\mathrm{keV}\ {\mathrm{cm}}^{2}$ (corresponding to ${t}_{\mathrm{cool}}\approx 200\,\mathrm{Myr}$ for the adopted cooling function) and $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$. Even though the central cooling time is $\ll 1\,\mathrm{Gyr}$, the system remains almost static for $\approx 2\,\mathrm{Gyr}$. Then a large condensation event raises the central entropy to $\sim 20\,\mathrm{keV}\ {\mathrm{cm}}^{2}$, which corresponds to $\min ({t}_{\mathrm{cool}})\approx 600\,\mathrm{Myr}$ and $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 15$. Afterward, the system remains relatively stable for ∼3 Gyr. The initial delay time before condensation ($\sim 10\,{t}_{\mathrm{cool}}$) is consistent with a period of buoyancy damping followed by steady exponential growth, as illustrated in Figure 3.

A similar interpretation suits run C1 from Sharma et al. (2012b), which starts with central entropy $\approx 1\,\mathrm{keV}\ {\mathrm{cm}}^{2}$ and $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 1$. Condensation begins promptly in that simulation and raises the ambient specific entropy at small radii to $\sim 10\,\mathrm{keV}\ {\mathrm{cm}}^{2}$ in $\lesssim 2\,\mathrm{Gyr}$, at which time $\min ({t}_{\mathrm{cool}})\,\sim 200\,\mathrm{Myr}$ and $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\sim 10$. Condensation then subsides, and the entropy and cooling-time profiles remain relatively stable for the next ∼2 Gyr, presumably for the same reason that run C10 resisted condensation for the first 2 Gyr.

In contrast, the transition to $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\gtrsim 10$ requires only ∼40 Myr for the Gaspari et al. (2013) simulations, which also enforce thermal balance within equipotential layers. There are three reasons for the action to proceed on a shorter timescale than in Sharma et al. (2012b). First, the greater metallicity (i.e., solar values) and lower temperature compared with Sharma et al. (2012b) result in a shorter central cooling time at a given central entropy level. Second, the gravitational potential is deeper because the gravity of the central galaxy's stars and a central supermassive black hole are both included, meaning that the freefall time at each radius is shorter. Third, and most importantly, Gaspari et al. (2013) continually drive turbulence, which counteracts the effects of buoyancy damping and raises the amplitude of the seed perturbations that thermal instability amplifies.

Taken together, the simulations of Sharma et al. (2012b) and Gaspari et al. (2013) have been viewed as support for the suggestion that systems in global thermal balance naturally stabilize near $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$, as observed in more realistic simulations of cluster cores regulated by precipitation-fed bipolar outflows (Gaspari et al. 2012; Li et al. 2015; Prasad et al. 2015; Meece et al. 2016; Yang & Reynolds 2016). However, both the analysis of Section 4 and the analytical toy model for configuration changes presented in the Appendix show that there is nothing unique about ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\approx 10$ from the point of view of local thermal stability. Instead, the Sharma et al. (2012b) simulations remain near $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$ for long time periods because increases in ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ stretch the time interval between condensation events to multiple gigayears, even though $\min ({t}_{\mathrm{cool}})$ remains $\lesssim 1\,\mathrm{Gyr}$.

Relaxing the requirement that each equipotential layer remain in thermal balance while retaining the global thermal balance constraint leads to a somewhat different interpretation for the minimum value of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$. When this restriction is relaxed, the entropy profile of a globally balanced system can change with time in response to the spatial distribution of heat input. In systems that can be divided into an inner isentropic zone and an outer power-law zone (see Figure 1 and the Appendix), the ambient ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ ratio is lowest at the boundary between those zones. Therefore, the minimum value of $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})$ can depend on how heat input is distributed with radius.

To see how the structure of a globally balanced system responds to the radial distribution of feedback heating, consider how heating affects the evolution of the entropy-profile slope. The Lagrangian derivative of ${\alpha }_{K}(r,t)$ with respect to time is

$\begin{eqnarray}\displaystyle \frac{d{\alpha }_{K}}{{dt}} & = & \displaystyle \frac{\partial }{\partial t}\left(\displaystyle \frac{\partial \mathrm{ln}K}{\partial \mathrm{ln}r}\right)+v\cdot \displaystyle \frac{\partial }{\partial r}\left(\displaystyle \frac{\partial \mathrm{ln}K}{\partial \mathrm{ln}r}\right)\\ & = & \displaystyle \frac{\partial }{\partial \mathrm{ln}r}\left(\displaystyle \frac{d\mathrm{ln}K}{{dt}}\right)+\left(\displaystyle \frac{\partial \mathrm{ln}K}{\partial \mathrm{ln}r}\right)\displaystyle \frac{\partial }{\partial \mathrm{ln}r}\left(\displaystyle \frac{v}{r}\right),\end{eqnarray} \tag{ 42 }$

where v is the radial velocity of the Lagrangian shell at r. That radial velocity depends on how quickly the volume bounded by the shell is expanding or contracting due to heating or cooling. It can be expressed in terms of the volume-weighted average of the change in specific volume ($1/\rho$) within r:

$\begin{eqnarray}\displaystyle \frac{v}{r}=\displaystyle \frac{d\mathrm{ln}r}{{dt}} & = & -\displaystyle \frac{1}{3}\left\langle \displaystyle \frac{d\mathrm{ln}\rho }{{dt}}\right\rangle ,\\ \left\langle \displaystyle \frac{d\mathrm{ln}\rho }{{dt}}\right\rangle & \equiv & \displaystyle \frac{3}{4\pi {r}^{3}}{\displaystyle \int }_{0}^{r}\left(\displaystyle \frac{d\mathrm{ln}\rho }{{dt}}\right)4\pi {r}^{2}\,{dr}.\end{eqnarray} \tag{ 43 }$

The Lagrangian derivative of ${\alpha }_{K}$ can therefore be written as

$\begin{eqnarray}\displaystyle \frac{d{\alpha }_{K}}{{dt}} & = & \displaystyle \frac{\partial }{\partial \mathrm{ln}r}\left(\displaystyle \frac{d\mathrm{ln}K}{{dt}}\right)\\ & & +{\alpha }_{K}\displaystyle \frac{d\mathrm{ln}\rho }{{dt}}-{\alpha }_{K}\left\langle \displaystyle \frac{d\mathrm{ln}\rho }{{dt}}\right\rangle .\end{eqnarray} \tag{ 44 }$

In this equation, the first term represents the changes in ${\alpha }_{K}$ that arise from the radial gradient of the net heating rate. The second term represents changes that arise from the compression or rarefaction of the local entropy gradient. The third term represents changes that result from the expansion or contraction in radius of the Lagrangian shell under consideration.

We are interested in knowing whether ${\alpha }_{K}$ increases or decreases with time in response to the distribution of heating and cooling with radius. If ${\alpha }_{K}=0$, then the answer is relatively simple. The entropy slope steepens in regions where the net heating rate increases with radius and declines in regions where it decreases with radius. However, the answer becomes more complicated when ${\alpha }_{K}\ne 0$ and the other two terms come into play. To better understand how the additional terms affect the evolution of the system as a whole, we explore three specific radial distributions of heating and cooling in the following three subsections, and in a fourth subsection (Section 5.3.4) we bring those explorations to bear on interpreting $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$.

5.3.1. Two Regions: Central Heating, Distributed Cooling

5.3.2. Two Regions: Distributed Heating, Central Cooling

The opposite case, in which heating exceeds cooling at large radii and cooling exceeds heating at small radii, evolves in the opposite sense. In such a system, the radial gradient of $d\mathrm{ln}K/{dt}$ is positive and tends to increase the entropy slope ${\alpha }_{K}$. However, increases in ${\alpha }_{K}$ make the second and third terms in Equation (44) more important. The effects of those other two terms are clearest if heating is insignificant. Then, cooling causes ρ to increase with time in the central cooling-dominated zone on a timescale $\sim {t}_{\mathrm{cool}}$. It also causes the third term to increase more than the second term, because the average cooling time within a sphere of radius r is generally shorter than the average cooling time at r in a system with a positive entropy gradient. As a result, the third term limits the entropy gradient to be no greater than the value at which the third term offsets the other two.

One can estimate the limiting value of ${\alpha }_{K}$ by assuming that the system remains approximately isothermal as it loses entropy through radiative cooling, which is usually valid for gas confined by a nearly isothermal gravitational potential, as long as ${t}_{\mathrm{cool}}\gg {t}_{\mathrm{ff}}$. With that assumption, we obtain

$\begin{eqnarray}&&\displaystyle \frac{d{\alpha }_{K}}{{dt}}\,\approx \,\displaystyle \frac{3{\alpha }_{K}}{{t}_{\mathrm{cool}}}\left(\displaystyle \frac{1-{\alpha }_{K}}{2-{\alpha }_{K}}\right)\end{eqnarray} \tag{ 45 }$

if the volume average in Equation (43) is performed assuming ${\alpha }_{K}$ is constant with radius. This estimate shows that cooling produces an inflow that approaches a limiting entropy slope of ${\alpha }_{K}\approx 1$, corresponding approximately to $\rho \propto {r}^{-3/2}$ in gas that remains nearly isothermal, and is in alignment with the self-similar cooling-flow solution derived by Bertschinger (1989). In a more realistic cluster-core potential, deviations from isothermality can slightly steepen the entropy slope of a steadily cooling system to ${\alpha }_{K}\approx 1.2$, but it always remains close to unity (e.g., Voit 2011). Formally, Equation (45) also suggests that a cooling region initialized without an entropy gradient remains isentropic, because the cooling time in the isothermal approximation we have adopted is then independent of radius. In that limit, the effects of temperature gradients and gradients in heat input can no longer be ignored. However, even if there is some growth of the entropy gradient, it proceeds on a timescale comparable to or longer than the cooling time.

We conclude that if cooling dominates heating at small radii, then entropy losses and the resulting compression cause the isentropic region to shrink in radius as the center of the system evolves toward a state with ${\alpha }_{K}\approx 1$. Condensation then happens only at small radii, where ${t}_{\mathrm{cool}}\lesssim {t}_{\mathrm{ff}}$. The result is to focus condensation and mass deposition onto the neighborhood of the black hole (e.g., Li & Bryan 2012). This outcome is highly favorable for the coupling of black-hole accretion to the global evolution of the surrounding medium. At larger radii, where heating exceeds cooling, the response of ${\alpha }_{K}$ depends on the radial gradient of heat input, which depends in turn on the feedback delivery mechanism. However, the resulting changes in ${\alpha }_{K}$ typically proceed on a much longer timescale than the cooling-induced changes to the inner profile, because the heat input required for global balance is spread over a much larger volume, in gas with much lower density and a much longer cooling time.

5.3.3. Three Regions: Heating–Cooling–Heating

Neither of the simple two-region cases yields a satisfactory explanation for how feedback in a cluster core can maintain a central cooling time $\lesssim 1\,\mathrm{Gyr}$ along with $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$ for periods of several gigayears. We must therefore take one more small step toward complexity and consider what happens if there are three nested regions in which the sign of $d\mathrm{ln}K/{dt}$ alternates, with heating exceeding cooling at the smallest and largest distances from the black hole, while cooling exceeds heating at intermediate radii.

Figure 4 schematically shows how the response to feedback of this three-region system contrasts with the two-region cases we have just discussed. The radial decline in $d\mathrm{ln}K/{dt}$ across the boundary between the innermost two regions causes convection that drives the inner part of the system toward an isentropic configuration that is prone to condensation. Meanwhile, the rise in $d\mathrm{ln}K/{dt}$ across across the boundary between the outermost two regions maintains a large-scale entropy gradient that suppresses thermal instability through buoyancy damping. The natural outcome of such a feedback response is a long-lasting configuration with an inner isentropic zone and an outer power-law zone (see Figure 1 again). In this configuration, the central production of cold clouds can fuel a feedback response that prevents condensation in the outer zone. The Appendix presents a toy model that makes this conceptual argument more quantitative.

Zoom In Zoom Out Reset image size

5.3.4. Interpretation of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\approx 10$

Figure 5 shows the feedback cycles in the Li et al. (2015) simulation. At the beginning, the gas is homogeneous at each radius and remains homogeneous as it cools, until condensation triggers strong feedback at $t=0.27\,\mathrm{Gyr}$. The amount of gas that is present within the central 20 kpc prior to the feedback outburst can be estimated by noting that ${\alpha }_{K}\approx 1$ and $K\approx 10\,\mathrm{keV}\ {\mathrm{cm}}^{2}$ at 20 kpc, giving a gas-mass estimate $\sim 3\times {10}^{11}\,{M}_{\odot }$ for a gas temperature ∼3 keV. Approximately $\sim 4\times {10}^{9}\,{M}_{\odot }$ of this gas condenses within ∼200 Myr and triggers feedback that leaves the system in a steady feedback-regulated state with $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 10$ for the next ∼1 Gyr, until star formation uses up all of the cold gas.

Zoom In Zoom Out Reset image size

Two more outburst and shutdown cycles follow the first one, and each cycle proceeds similarly. Without heat input, the core cools homogeneously and develops a power-law central entropy gradient that focuses condensation on the black hole. When the black hole erupts, uplift stimulates rapid condensation and spreads condensing gas over a much larger volume. After each eruption, the system is left in a long-lasting, well-regulated state with a flat central entropy gradient.

In order to understand why these cycles unfold as they do, it is important to look closely at the distribution of specific entropy K at each radius. Figure 6 shows the distribution of K as a function of r at t = 1 Gyr, which is in the midst of the first feedback cycle. Broadly speaking, the entropy profile is approximately isentropic in the inner regions and approximately a power law with ${\alpha }_{K}\approx 1$ in the outer regions. The transition from ${\alpha }_{K}\ll 1$ to ${\alpha }_{K}\approx 1$ is at the radius where ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\approx 10$. However, gas at any given radius can span a wide range in K.

Zoom In Zoom Out Reset image size

The actual dispersion in specific entropy is not quite as broad as the figure's logarithmic heat map makes it seem. Lines in the inset figure show the median entropy at each radius (solid line) and the 20%–80% percentile range (dashed lines). Despite the wide range of extremes in entropy, the dispersion is actually quite narrow at large radii. The extremes are associated with the high-entropy outflow and the lower-entropy gas uplifted by it (Li & Bryan 2014b), but buoyancy damping keeps most of the gas nearly homogeneous at radii where the entropy gradient is significant.

In the center, where the entropy gradient vanishes, buoyancy damping is less successful at suppressing entropy perturbations, and thermal instability is able to produce condensation (Section 3). Gas blobs below the median entropy in this zone are declining in specific entropy and will soon join the condensed component. Outside of this thermally unstable zone, condensation is scarcer. At the epoch shown, there are two isolated patches of multiphase gas at larger radii, roughly 30 and 50 kpc from the center. Buoyancy damping cannot suppress condensation in these patches because they have been uplifted by the outflow (Section 2) and are on nearly ballistic trajectories.

The six panels of Figure 7 show how the entropy distribution evolves during the first outburst of feedback. Before the outburst erupts at $t=0.27\,\mathrm{Gyr}$, the K distribution at each r is still nearly homogeneous for two reasons: (1) pure cooling without compensating heating does not lead to phase separation (Section 5.1), and (2) even if there were some entropy perturbations, they would be damped by buoyancy (Section 3.4). The next panel (t = 0.37 Gyr) shows the system 100 Myr after the beginning of the outburst, and the leading edge of the high-entropy outflow (at ∼40 kpc) is well defined. Gas at larger radii is still undisturbed, and the median entropy still follows an approximate power law to $\ll 10\,\mathrm{kpc}$. However, there are significant amounts of condensing low-entropy gas at ∼20 kpc. That cannot be gas condensing out of the ambient medium at 20 kpc, because the median cooling time there was $\gtrsim 200\,\mathrm{Myr}$ at the time the outburst began. Instead, it is low-entropy gas that was transported outward from smaller radii by the outflow, and uplift is allowing it to condense.

Zoom In Zoom Out Reset image size

After another 100 Myr (at t = 0.47 Gyr), an isentropic central region has started to develop at ${K}_{0}\approx 10\,\mathrm{keV}\ {\mathrm{cm}}^{2}$, where the cooling time is $\approx 200\,\mathrm{Myr}$. In order for an isentropic central region to form, the bulk of the lower-entropy gas must be eliminated through either condensation, heating, or mixing. Figure 5 shows that less than $\sim {10}^{10}\,{M}_{\odot }$ condenses during this event, implying that the rest of the $\sim 3\times {10}^{11}\,{M}_{\odot }$ initially within the central $\lesssim 20\,\mathrm{kpc}$ is either heated or mixed into gas of greater entropy.

The remaining three frames, spanning the next 300 Myr, show that the system becomes well-regulated after arriving at a state with an isentropic core and a power-law outer region (see Section 5.3.3). Condensation continues in the isentropic core at $\lesssim 20\,\mathrm{kpc}$ but proceeds only in isolated patches beyond where the entropy gradient steepens, and those patches are moving on approximately ballistic trajectories. The outflow continues to populate the extremes of the entropy distribution. However, buoyancy damping causes those extremes to evolve toward the median, as shown by the lines in the lower-right panel, which mark the median and 20%–80% percentile range, as in the Figure 6 inset.

Figure 8 shows six more frames, this time spanning the period from the end of the first feedback episode to a time approximately 400 Myr after the beginning of the second. At $t=1.44\,\mathrm{Gyr}$, the cold fuel has just run out, causing the AGN to shut down, but the residual outflow has not finished propagating. Some of the gas is far from the median entropy level, and the central cooling time is $\approx 200\,\mathrm{Myr}$. During the next ∼100 Myr, buoyancy damping reduces the extremes in entropy at each radius in the power-law zone. Meanwhile, the central cooling time begins to decline, but it does so almost uniformly, because there is no source of free energy to cause phase separation (Section 5.1). Condensation finally starts to happen as the simulation approaches $t=1.59\,\mathrm{Gyr}$, and it develops most quickly at the boundary between the isentropic and power-law zones, because that is where ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ is smallest (Section 5.3).

Zoom In Zoom Out Reset image size

Once condensation begins, cooling gas quickly converges toward the center and triggers a second feedback outburst at $t=1.64\,\mathrm{Gyr}$. At this moment, the median entropy distribution is nearly a power law with ${\alpha }_{K}\approx 1.3$ from 100 kpc to within 10 kpc. This median slope is almost identical to the triggering slope at $t=0.37\,\mathrm{Gyr}$ in Figure 7, but the pre-triggering entropy distribution at $\lesssim 10\,\mathrm{kpc}$ is much less homogeneous because of the condensation that started happening 50 Myr earlier.

A major consequence of this pre-existing inhomogeneity is that uplift induces condensation far more effectively when the outflow begins, yielding $\sim 5\times {10}^{10}\,{M}_{\odot }$ of cold gas during the next 200 Myr (see Figure 5). The gas mass within ∼20 kpc of the center as condensation begins is similar to the $\sim 3\times {10}^{11}\,{M}_{\odot }$ that was there during the first outburst, but 10 times more of it condenses (see Figure 5). The rest must be either heated or mixed to $\gtrsim 20\,\mathrm{keV}\ {\mathrm{cm}}^{2}$ because the central region becomes nearly isentropic at this level by $t=2.05\,\mathrm{Gyr}$, as shown in the final panel of Figure 8.

After the second outburst produces another reservoir of cold gas fuel, the system enters a long period of steady self-regulation. The median entropy profile within ∼20 kpc stays flat at ${K}_{0}\approx 20\,\mathrm{keV}\ {\mathrm{cm}}^{2}$, so that $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 15$, apart from some stochastic excursions as the black-hole energy output fluctuates. However, the cold fuel eventually runs out after sustaining $\approx 2.5\,\mathrm{Gyr}$ of activity, setting the stage for the third outburst.

Figure 9 concentrates on the events leading up to the third outburst, which require more time to unfold because the central cooling time of the system is $\approx 600\,\mathrm{Myr}$ when the cold fuel runs out. Complete depletion of cold fuel could be a consequence of the idealizations in this simulation that never happen in real cluster cores. We mention this possibility because very few real clusters have central cooling times $\lesssim 1\,\mathrm{Gyr}$ and no multiphase gas (Voit & Donahue 2015). However, the period of quiescence before the third outburst in this simulation beautifully illustrates the physics of AGN triggering via condensation, and similar processes are likely to operate in less idealized systems.

Zoom In Zoom Out Reset image size

The events that lead to triggering are almost identical to those in Figure 8, but here the progression of thermal instability is easier to see. Buoyancy damping again reduces entropy perturbations in the power-law zone, and the central regions again cool uniformly because there is no source of free energy to promote phase separation. However, buoyancy damping is least effective at the junction between the isentropic and power-law zones, where ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ is smallest and is also declining with time. That junction moves inward as the isentropic zone loses thermal energy and becomes more highly compressed. As a result, the median entropy profile approaches a pure power law in radius, while thermal instability progresses toward condensation at the outer edge of the isentropic zone. Condensation then occurs between $t=4.77\,\mathrm{Gyr}$ and $t=4.82\,\mathrm{Gyr}$. It happens first at $r\approx 20\,\mathrm{kpc}$, but soon the entire region within 20 kpc contains multiphase gas.

The feedback outflow then sweeps through the condensing region, and 200 Myr later (t = 5.00 Gyr) a new isentropic zone has begun to form. This triggering period produces even more condensation than the previous two outbursts, yielding $\approx {10}^{11}\,{M}_{\odot }$ of cold gas. With this fuel source, the feedback outflow stabilizes the system's global configuration at ${K}_{0}\approx 30\,\mathrm{keV}\ {\mathrm{cm}}^{2}$ and $\min ({t}_{\mathrm{cool}}/{t}_{\mathrm{ff}})\approx 20$ for another ∼1 Gyr.

During the long periods of relatively steady feedback regulation between these triggering events, condensation happens mostly within the isentropic zone, where it is not inhibited by buoyancy damping. However, multiphase regions are sometimes observed far outside the isentropic zone, at radii where ${\alpha }_{K}\approx 1$ and ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\gt 10$ in the ambient medium, implying that buoyancy damping should be suppressing condensation. In plots of the K(r) distribution, those multiphase regions show up as isolated vertical stripes extending far below the median entropy level.

The origin of this multiphase gas in the power-law zone is not the growth of linear thermal instability but rather the uplift. These gas patches originate in the isentropic zone, are uplifted by the outflow, and then follow trajectories that are approximately ballistic. Black arrows in Figure 10 show how one such region rises and falls. It can be seen rising out of the isentropic zone, reaching $\approx 30\,\mathrm{kpc}$ at $t=5.10\,\mathrm{Gyr}$. Between $t=5.20\,\mathrm{Gyr}$ and $t=5.30\,\mathrm{Gyr}$, it reaches an apex at $\approx 80$–90 kpc and falls back through $\approx 30\,\mathrm{kpc}$ at $t\approx 5.42\,\mathrm{Gyr}$. If the motion of the gas blob were purely ballistic, this excursion in radius would take $\approx 240\,\mathrm{Myr}$ and would require the projectile to rise and fall through r = 30 kpc, where the circular velocity is $\approx 550\,\mathrm{km}\ {{\rm{s}}}^{-1}$, at a speed $\gtrsim 900\,\mathrm{km}\ {{\rm{s}}}^{-1}$. Instead, the peak speed during both the upward and downward portions of the trajectory is $\sim 500\,\mathrm{km}\ {{\rm{s}}}^{-1}$, while most of the time is spent at smaller speeds near the trajectory's turning point. These deviations from ballistic motion imply that the outflow is still providing impetus to the multiphase gas blob during the upward stage, and that drag is limiting the infall speed to $\approx 550\,\mathrm{km}\ {{\rm{s}}}^{-1}$ during the downward stage.

Zoom In Zoom Out Reset image size

It can be argued that these multiphase gas clumps would be deviating more strongly from ballistic trajectories in a more realistic and highly resolved simulation because of hydrodynamic shredding and drag (e.g., Nulsen 1986; Pizzolato & Soker 2005). In this simulation, however, drag is clearly not necessary to promote the development of multiphase gas at large altitudes, because the gas is multiphase on both the upward and downward portions of the trajectory. High-resolution simulations specifically designed to track the hydrodynamical evolution (and perhaps even magnetohydrodynamical evolution) of uplifted gas will be needed to draw more certain conclusions about the roles of uplift and drag in promoting condensation and multiphase structure within the power-law zone.

The rotation of the system, which we have ignored, has at least two major consequences: (1) locally it can reduce the effective gravitational potential and therefore the stratification on which buoyancy damping depends, and (2) globally it can inhibit the accretion of condensed gas onto the central black hole. Consequence (1) can promote condensation in regions where it otherwise would not happen. Consequence (2) has the potential to decouple condensation from black-hole fueling. However, Gaspari et al. (2015) have explored how rotation affects condensation and chaotic cold accretion using high-resolution simulations of thermally balanced media and show that precipitation can still deliver plenty of cold fuel to the vicinity of the black hole. The key is for turbulence and cloud−cloud collisions to continually populate the low end of the distribution of specific angular momentum among cold clouds. This is more easily done when the velocity dispersion of turbulence exceeds the mean rotational speed.

It is tempting to speculate about the implications of consequence (1) for the replenishment of condensed gas in galactic disks. In a rotating frame of reference, the effective gravitational potential of a spherically symmetric halo with a stellar disk becomes more like that in the plane-parallel simulations of Meece et al. (2015). If the halo contains hydrostatic ambient gas with ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\gg 1$, then buoyancy damping can inhibit condensation everywhere except the midplane. Condensation of hot ambient halo gas may then replenish the disk's reservoir of cold gas without any obvious evidence for infall in the form of cool clouds at larger altitudes. For a more formal analysis of thermal instability in a rotating circumgalactic medium, see Nipoti & Posti (2014).

If cooling of the circumgalactic medium depends strongly on metallicity, as it certainly will in systems with virial temperatures $\lesssim 1\,\mathrm{keV}$, then inhomogeneities in the distribution of metals will promote thermal instability and condensation of the more highly enriched gas (e.g., Fraternali & Binney 2008; Marinacci et al. 2010; Oppenheimer et al. 2010). An over-enriched gas blob with the same K and P as the medium surrounding it will cool faster, causing it to sink, and it will not respond as strongly to buoyancy damping. This coupling among enrichment, entropy, and buoyancy has complex consequences for condensation, meaning that metal transport and mixing in circumgalactic gas will be important to simulate properly in order to clarify the role of inhomogeneous enrichment.

Perhaps counterintuitively, inhomogeneities in enrichment could make it easier for lower-metallicity gas to condense within the galactic disk itself. This may happen if high-metallicity gas more easily condenses into cold, ballistic clouds at greater altitudes, while lower metallicity gas settles more homogeneously into the disk and condenses there instead. Again, high-resolution targeted simulations including metal transport will be necessary to resolve this question.

Our entire analysis has implicitly assumed the existence of tangled magnetic fields because we have ignored electron thermal conduction (see Soker 2010b for a discussion of the role of magnetic fields in cooling flows). In reality, conduction is likely to be somewhat anisotropic without being completely suppressed, which will naturally cause condensation to be filamentary (e.g., Sharma et al. 2010; McCourt et al. 2012; Wagh et al. 2014). Also, larger-scale conduction may help lift cluster cores out of precipitation-regulated states, if the central entropy rises enough for radial heat transport to exceed radiative cooling in the core (Guo et al. 2008; Voit et al. 2008, 2015b; Parrish et al. 2010). If that happens, then feedback is no longer necessary to prevent catastrophic cooling.

Another potential consequence of magnetic fields is enhanced drag, which can promote condensation (Section 3.3.1). This has long been recognized as a possibility (e.g., Nulsen 1986), and recent ALMA observations of remarkably small, sub-Keplerian velocity dispersions among the bulk of the molecular gas clouds in galaxy-cluster cores (McNamara et al. 2014; Russell et al. 2016) have rekindled interest in it (McNamara et al. 2016). Obviously, magnetohydrodynamical simulations are required to clarify whether magnetic tension is strong enough to suspend molecular clouds within a medium that is many orders of magnitude less dense.

Cosmological MHD simulations will also be needed to clarify whether the magnetothermal instability (Balbus 1991) or heat-flux driven buoyancy instability (Quataert 2008; Parrish & Quataert 2008) alter the general picture of susceptibility to condensation that we have outlined here. When anisotropic conduction is important, buoyancy is more closely related to temperature gradients than to entropy gradients (e.g., Balbus 2000; Sharma et al. 2009), which is likely to reduce the dependence of buoyancy damping on ${\alpha }_{K}$. Recent simulations have shown that maximally efficient anisotropic conduction does indeed alter the effects of AGN feedback on the intracluster medium (Kannan et al. 2017).

We have emphasized the importance of heat transport in stabilizing the circumgalactic medium of massive galaxies but have not paid any attention to how the transported energy is transformed into heat. The dissipation of feedback energy in galaxy-cluster cores is currently an active area of research. Many possibilities are being investigated, including turbulence, shock heating, and mixing of hot gas with cooler gas (e.g., Banerjee & Sharma 2014; Zhuravleva et al. 2014; Hillel & Soker 2016; Yang & Reynolds 2016). These dissipation processes will need to be included in a more complete model.

The temperature and electron density profiles in the isentropic region are given by

$\begin{eqnarray}T(r) & = & {T}_{0}\left[1-\displaystyle \frac{3{\alpha }_{K}}{5}\mathrm{ln}\left(\displaystyle \frac{r}{{r}_{0}}\right)\right],\\ {n}_{e}(r) & = & {\left(\displaystyle \frac{{{kT}}_{0}}{{K}_{0}}\right)}^{3/2}{\left[1-\displaystyle \frac{3{\alpha }_{K}}{5}\mathrm{ln}\left(\displaystyle \frac{r}{{r}_{0}}\right)\right]}^{3/2},\end{eqnarray} \tag{ 48 }$

and integrating to obtain the ambient gas mass inside r₀ yields

$\begin{eqnarray}{M}_{0} & = & \displaystyle \frac{4\pi }{3}{r}_{0}^{3}\,{\mu }_{e}{m}_{p}{\left(\displaystyle \frac{4{{kT}}_{\phi }}{3{K}_{0}}\right)}^{3/2}\displaystyle \frac{{I}_{3/2}({\alpha }_{K})}{{\alpha }_{K}^{3/2}},\\ {I}_{n}({\alpha }_{K}) & \equiv & {\displaystyle \int }_{0}^{1}{\left(1-\displaystyle \frac{3{\alpha }_{K}}{5}\mathrm{ln}x\right)}^{n}3{x}^{2}{dx},\end{eqnarray} \tag{ 49 }$

with ${I}_{3/2}({\alpha }_{K})\approx (1+0.33\,{\alpha }_{K}^{1.1})$ to within 0.2% in the range $2/3\lt {\alpha }_{K}\lt 2$. Reorganizing the equation for M₀ then gives r₀ in terms of the three adjustable parameters (K₀, M₀, ${\alpha }_{K}$) that control the system's configuration.

In most astrophysical systems of interest, the ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ ratio declines with radius in isentropic gas. In this particular system, the dependence of cooling time on radius in the isentropic region follows

$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}{t}_{\mathrm{cool}}}{d\mathrm{ln}r}=(1+2\lambda )\displaystyle \frac{3{\alpha }_{K}}{10}{\left[1-\displaystyle \frac{3{\alpha }_{K}}{5}\mathrm{ln}\left(\displaystyle \frac{r}{{r}_{0}}\right)\right]}^{-1},\end{eqnarray} \tag{ 50 }$

where $\lambda \equiv d\mathrm{ln}\,{\rm{\Lambda }}/d\mathrm{ln}T$. The condition for ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ to decline with radius in a system with constant T_ϕ is $d\mathrm{ln}{t}_{\mathrm{cool}}/d\mathrm{ln}r\lt 1$. A decline with radius is therefore guaranteed for ${\alpha }_{K}\lt 5/3$, since astrophysical cooling functions have $\lambda \leqslant 1/2$ at the temperatures characteristic of ambient circumgalactic gas.

Buoyancy damping does not suppress thermal instability in this part of the system. It is therefore prone to producing multiphase gas via condensation on a timescale $\sim \mathrm{ln}(1/{\delta }_{K})\times {t}_{\mathrm{cool}}$, where ${\delta }_{K}$ is the fractional amplitude of the initial entropy fluctuations.

The structure of the system's isothermal region outside of r₀ is given by

$\begin{eqnarray}K(r) & = & {K}_{0}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{{\alpha }_{K}},\\ {n}_{e}(r) & = & {\left(\displaystyle \frac{{{kT}}_{0}}{{K}_{0}}\right)}^{3/2}{\left(\displaystyle \frac{r}{{r}_{0}}\right)}^{-3{\alpha }_{K}/2},\\ \displaystyle \frac{{t}_{\mathrm{cool}}}{{t}_{\mathrm{ff}}} & \propto & {r}^{(3{\alpha }_{K}-2)/2}.\end{eqnarray} \tag{ 51 }$

Buoyancy damping suppresses condensation in the isothermal region as long as ${\alpha }_{K}\gtrsim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$. Suppression of condensation therefore applies to the entire isothermal region if ${\alpha }_{K}\gt 2/3$ and ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}\gg 1$ at r₀. In such systems, the minimum value of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ is at r₀.

In order for the system to remain in overall thermal balance as it evolves, the heat energy lost by one of its subsystems must be gained by another subsystem. Let subsystem 1 with gas mass $\delta {M}_{g}$ be the one that loses an amount of heat energy $\delta Q$ and experiences a corresponding change in specific entropy from K₁ to ${K}_{1}-\delta {K}_{1}$, such that

$\begin{eqnarray}&&\delta Q=\displaystyle \frac{3k}{2}\displaystyle \frac{\delta {M}_{g}}{{\mu }_{e}{m}_{p}}{\int }_{{K}_{1}-\delta {K}_{1}}^{{K}_{1}}{T}_{1}(K)\displaystyle \frac{{dK}}{K}.\end{eqnarray} \tag{ 52 }$

In general, heat transfer changes the temperature T₁ of the subsystem, and the relationship between $\delta Q$ and $\delta K$ depends on ${T}_{1}(K)$ during the heat-transfer episode. For example, the particular case of condensation at constant pressure corresponds to ${T}_{1}\propto {K}^{\eta }$ with $\eta =3/5$, and condensation at constant volume corresponds to $\eta =1$. In either case, a condensation process in which ${K}_{1}-\delta {K}_{1}\ll {K}_{1}$ leads to

$\begin{eqnarray}&&\delta Q=\displaystyle \frac{3}{2}\displaystyle \frac{{{kT}}_{1}(K)}{\eta }\displaystyle \frac{\delta {M}_{g}}{{\mu }_{e}{m}_{p}}.\end{eqnarray} \tag{ 53 }$

A more complex condensation process can be characterized by the effective value of η for which this equation holds.

Let subsystem 2 be that one that gains heat energy $\delta Q$. If subsystem 2 has a gas mass ${M}_{g}\gg \delta {M}_{g}$, then the gain in heat energy raises its specific entropy by

$\begin{eqnarray}&&\delta {K}_{2}=\displaystyle \frac{{K}_{2}}{\eta }\displaystyle \frac{{T}_{1}({K}_{1})}{{T}_{2}({K}_{2})}\displaystyle \frac{\delta {M}_{g}}{{M}_{g}}.\end{eqnarray} \tag{ 54 }$

Note that this equation can also be used to derive the relation ${M}_{g}\propto {K}^{-\eta }$ between the specific entropy K and gas mass M_g of the ambient medium in a system that is kept in global thermal balance as condensation gradually removes mass from the ambient phase and adds it to the condensed phase.

Heat redistribution algorithms that keep each equipotential layer in independent thermal balance (e.g., McCourt et al. 2012; Sharma et al. 2012b; Gaspari et al. 2013; Meece et al. 2015) allow gas to condense out of the ambient phase in the isentropic region and produce an accompanying rise in K₀. Meanwhile, the specific entropy K of each Lagrangian layer outside of the isentropic region remains unchanged because buoyancy damping causes entropy fluctuations there to saturate at a fractional amplitude ${\delta }_{K}\sim {({t}_{\mathrm{ff}}/{t}_{\mathrm{cool}})}^{2}$ instead of condensing. Consequently, a rise in K₀ drives convective mixing in the neighborhood of r₀ as gas with $K\lt {K}_{0}$ sinks into the isentropic zone and adds to M₀.

Those changes in system configuration are best analyzed in the gas-mass domain. Let M_g (K) represent the mass of ambient gas with specific entropy less than K, and let M_c be the mass of gas that has already condensed. Inverting M_g (K) gives a two-segment entropy distribution that is isentropic at K₀ for ${M}_{g}\lt {M}_{0}$ and follows a power law with $d\mathrm{ln}K/d\mathrm{ln}{M}_{g}=2{\alpha }_{K}/(6-3{\alpha }_{K})$ for ${M}_{g}\gt {M}_{0}$. A condensation event that raises K₀ by $\delta {K}_{0}$ adds gas mass $\delta {M}_{c}=\eta {M}_{0}\,\delta {K}_{0}/{K}_{0}$ to the condensed phase and changes M₀ by

$\begin{eqnarray}&&\delta {M}_{0}\ =\ {M}_{0}\left(\displaystyle \frac{6-3{\alpha }_{K}}{2{\alpha }_{K}}-\eta \right)\displaystyle \frac{\delta {K}_{0}}{{K}_{0}}.\end{eqnarray} \tag{ 55 }$

Condensation therefore drives the ratio of condensed to isentropic gas toward ${M}_{c}/{M}_{0}=2{\alpha }_{K}\eta /(6-3{\alpha }_{K})$, implying that the condensed gas mass is comparable to M₀ for representative values of ${\alpha }_{K}\sim 1$ and $\eta \approx 3/5$.

As K₀ rises, there is a particular value of η that ensures that the K(r) profile remains isothermal,

$\begin{eqnarray}{\eta }_{0} & = & \displaystyle \frac{3}{2}-\displaystyle \frac{3}{{\alpha }_{K}}\left[1-\displaystyle \frac{1}{{I}_{3/2}({\alpha }_{K})}\right]\\ & \approx & \displaystyle \frac{3}{2}\left[\displaystyle \frac{3+{\alpha }_{K}^{1.1}-2{\alpha }_{K}^{0.1}}{3+{\alpha }_{K}^{1.1}}\right],\end{eqnarray} \tag{ 56 }$

because it ensures that ${M}_{0}+{M}_{c}={M}_{0}[2/(2-{\alpha }_{K})]{I}_{3/2}^{-1}({\alpha }_{K})$, which is the gas mass that would be contained within r if the isothermal region extended all the way to r = 0. For ${\alpha }_{K}\approx 1$, this equation gives ${\eta }_{0}\approx 3/4$, which is in between the characteristic values for constant-pressure condensation ($\eta =3/5$) and constant-volume condensation ($\eta =1$). Evolution of K(r) in a system with thermally balanced equipotential layers therefore remains close to the idealized isentropic/isothermal configuration, with K(r) staying nearly fixed at $r\gt {r}_{0}$ as condensation in lower-lying layers raises the value of K₀ (see, e.g., Sharma et al. 2012b; Gaspari et al. 2013). Deviations of the actual value of η from ${\eta }_{0}$ cause M₀ to deviate somewhat from the required value for maintaining the idealized configuration and result in small deviations from isothermality just outside of r₀.

The rate at which such a system's K(r) profile evolves depends primarily on ${t}_{\mathrm{cool}}({K}_{0})$ but also on the initial sizes of the entropy fluctuations that lead to condensation. If the fractional amplitudes of these fluctuations are of order unity, then $d\mathrm{ln}{K}_{0}/{dt}\approx {t}_{\mathrm{cool}}^{-1}({K}_{0})$. If instead the seed fluctuations are smaller, then more time is required for condensation. The results of Section 4 show that buoyancy damping leads to saturation of perturbation growth at smaller amplitudes for larger values of ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$. Condensation at the center of the system, where the entropy gradient vanishes, is still possible but requires a time period that grows progressively longer compared to the cooling time as ${t}_{\mathrm{cool}}/{t}_{\mathrm{ff}}$ increases.

A globally balanced system in which one layer gains heat energy while another layer loses the same amount of heat energy has more freedom to change configuration. Qualitatively, the two simplest kinds of configuration changes are (1) if cooling exceeds heating inside of r₀, causing K₀ to decrease, then heating must exceed cooling in the outer region, or (2) if heating exceeds cooling inside of r₀, causing K₀ to increase, then cooling must exceed heating in the outer region.

In our toy model, a small change in entropy $\delta K({M}_{g})$ as a function of the gas mass M_g enclosed within a given radius corresponds to a change in heat energy given by

$\begin{eqnarray}\delta Q & = & \displaystyle \frac{3}{2}\displaystyle \frac{{{kT}}_{0}}{\mu {m}_{p}}\left[\delta \mathrm{ln}{K}_{0}\cdot {\displaystyle \int }_{0}^{{M}_{0}}\displaystyle \frac{T}{{T}_{0}}\,{{dM}}_{g}\right.\\ & & \left.+{\displaystyle \int }_{{M}_{0}}^{{M}_{{\rm{b}}}}\delta \mathrm{ln}K\cdot {{dM}}_{g}\right],\end{eqnarray} \tag{ 57 }$

where ${M}_{{\rm{b}}}$ is the ambient gas mass within the boundary radius ${r}_{{\rm{b}}}$. The configuration changes permitted by global thermal balance within the radius ${r}_{{\rm{b}}}$ have the property

$\begin{eqnarray}&&{\displaystyle \int }_{{M}_{0}}^{{M}_{{\rm{b}}}}\delta \mathrm{ln}K({M}_{g})\cdot {{dM}}_{g}\\ &&\quad =\,-{M}_{0}\,\displaystyle \frac{{I}_{5/2}({\alpha }_{K})}{{I}_{3/2}({\alpha }_{K})}\,\delta \mathrm{ln}{K}_{0},\end{eqnarray} \tag{ 58 }$

where ${I}_{5/2}({\alpha }_{K})$ and ${I}_{3/2}({\alpha }_{K})$ are both given by the definition in Equation (49). From Equations (49) and (51), we obtain an expression for the entropy-profile change in terms of changes in the three model parameters:

$\begin{eqnarray}\delta \mathrm{ln}K & = & \left(1-\displaystyle \frac{{\alpha }_{K}}{2}\right)\delta \mathrm{ln}{K}_{0}-\displaystyle \frac{{\alpha }_{K}}{3}\,\delta \mathrm{ln}{M}_{0}\\ & & -\left[\mathrm{ln}\left(\displaystyle \frac{K}{{K}_{0}}\right)-\displaystyle \frac{{\alpha }_{K}}{2}+\displaystyle \frac{{\alpha }_{K}}{3}\displaystyle \frac{d\mathrm{ln}{I}_{3/2}}{d\mathrm{ln}{\alpha }_{K}}\right]\delta \mathrm{ln}{\alpha }_{K}.\end{eqnarray} \tag{ 59 }$

Integrating this function over M_g then gives

$\begin{eqnarray}&&\left(\displaystyle \frac{{M}_{b}-{M}_{0}}{{M}_{0}}\right)\,\delta \mathrm{ln}{K}_{b}-\left(\displaystyle \frac{2{\alpha }_{K}}{6-3{\alpha }_{K}}\right)\left[\left(\displaystyle \frac{{M}_{b}-{M}_{0}}{{M}_{0}}\right)\right.\\ &&\quad \left.+\displaystyle \frac{3}{{\alpha }_{K}{I}_{3/2}}\right]\delta \mathrm{ln}{\alpha }_{K}=-\displaystyle \frac{{I}_{5/2}}{{I}_{3/2}}\delta \mathrm{ln}{K}_{0},\end{eqnarray} \tag{ 60 }$

where $\delta \mathrm{ln}{K}_{b}$ is the change in specific entropy at the boundary, obtained by substituting ${K}_{b}\equiv K({r}_{{\rm{b}}})$ into Equation (59). Parameter adjustments satisfying this equation change the system's configuration without changing the amount of heat energy it contains.

Some configuration changes require a flow of free energy through the system, while others do not. The change in the system's global entropy S when K₀ changes under conditions of global thermal balance is

$\begin{eqnarray}\delta S & = & \displaystyle \frac{3k}{2\mu {m}_{p}}{\displaystyle \int }_{0}^{{M}_{{\rm{b}}}}\displaystyle \frac{\delta K}{K}\ {{dM}}_{g}\\ & = & -\displaystyle \frac{3{{kM}}_{0}}{2\mu {m}_{p}}\left[\displaystyle \frac{{I}_{5/2}}{{I}_{3/2}}-1\right]\delta \mathrm{ln}{K}_{0}.\end{eqnarray} \tag{ 61 }$

A reduction of K₀ can proceed without a source of free energy (since ${I}_{5/2}\gt {I}_{3/2}$) because it corresponds to a transfer of heat from the higher temperature isentropic zone to the lower temperature isothermal zone, thereby boosting the global entropy. However, a growth of K₀ requires a source of free energy, because the lower-temperature outer regions then lose heat while the higher-temperature inner regions gain it, resulting in $\delta S\lt 0$.

With the global balance constraint in place, the system has only two degrees of freedom. Applying another constraint at the boundary leaves only one degree of freedom. One option is to constrain the temperature at the outer boundary to remain at T₀, which has the consequence of holding ${\alpha }_{K}$ fixed. In that case, Equation (58) reduces to

$\begin{eqnarray}\delta \mathrm{ln}{M}_{0} & = & \left[\left(\displaystyle \frac{6-3{\alpha }_{K}}{2{\alpha }_{K}}\right)+\displaystyle \frac{3}{{\alpha }_{K}}\displaystyle \frac{{I}_{5/2}}{{I}_{3/2}}\right.\\ & & \left.\times \left(\displaystyle \frac{{M}_{0}}{{M}_{b}-{M}_{0}}\right)\right]\delta \mathrm{ln}{K}_{0},\end{eqnarray} \tag{ 62 }$

and central heating with ${\alpha }_{K}\lt 2$ causes both K₀ and M₀ to rise while ${\alpha }_{K}$ stays fixed. Alternatively, holding the entropy constant at the boundary gives

$\begin{eqnarray}\delta \mathrm{ln}{K}_{0} & = & \left(\displaystyle \frac{2}{2-{\alpha }_{K}}\right)\left[\displaystyle \frac{{I}_{3/2}}{{I}_{5/2}}\right.\\ & & \left.\times \left(\displaystyle \frac{{M}_{b}-{M}_{0}}{{M}_{0}}\right)+\displaystyle \frac{3}{{I}_{5/2}{\alpha }_{K}}\right]\delta \mathrm{ln}{\alpha }_{K},\end{eqnarray} \tag{ 63 }$

$\begin{eqnarray}\delta \mathrm{ln}{M}_{0} & = & \displaystyle \frac{3}{{\alpha }_{K}}\left[\displaystyle \frac{{I}_{3/2}}{{I}_{5/2}}\left(\displaystyle \frac{{M}_{b}-{M}_{0}}{{M}_{0}}\right)+\displaystyle \frac{3}{{I}_{5/2}{\alpha }_{K}}\right.\\ & & \left.-\mathrm{ln}\displaystyle \frac{{K}_{b}}{{K}_{0}}+\displaystyle \frac{{\alpha }_{K}}{2}-\displaystyle \frac{{\alpha }_{K}}{3}\displaystyle \frac{d\mathrm{ln}{I}_{3/2}}{d\mathrm{ln}{\alpha }_{K}}\right]\delta \mathrm{ln}{\alpha }_{K}.\end{eqnarray} \tag{ 64 }$

For this boundary condition, M₀ must again increase as K₀ grows in order for the system to remain in global thermal balance. The entropy slope ${\alpha }_{K}$ must increase as well, but that change is small if ${M}_{b}\gg {M}_{0}$.

This result has important qualitative implications for feedback-regulated systems. If they are thermally balanced by a central heat source, then central heating makes ambient gas there isentropic, and the values of both K₀ and M₀ increase with time. Condensation in that region is not suppressed by buoyancy damping and can add to the mass M_c of condensed gas at a rate ${\dot{M}}_{c}\sim {M}_{0}/{t}_{\mathrm{cool}}({K}_{0})$ until the rise in K₀ is sufficient for ${t}_{\mathrm{cool}}({K}_{0})$ to approach a Hubble time. This is not a finely regulated feedback loop but rather a monotonically evolving system in which feedback gradually shuts off condensation while a potentially large amount of condensed gas accumulates.

In contrast, globally balanced systems in which heating exceeds cooling outside of the isentropic zone generally have values of K₀ and M₀ that decrease with time. Buoyancy damping limits condensation to a relatively small central region, and the associated cooling time also decreases with time, with ${t}_{\mathrm{cool}}\mathop{\sim }\limits^{\propto }{K}_{0}^{3/2}$. Ultimately, something must prevent the central cooling time in the ambient medium from becoming singular. One lower limit on the cooling time of ambient gas at small radii is ${t}_{\mathrm{cool}}\gtrsim {t}_{\mathrm{ff}}$, because otherwise thermal instability is unavoidable. Another is stochastic heating that occasionally exceeds cooling, even if cooling exceeds heating during most of the duty cycle. In a less idealized system, the ambient cooling time that prevails after a heating outburst will be the timescale on which coupling and self-regulation happen, and this timescale can be much less than a Hubble time.

Stable global regulation at a particular K₀ and r₀ requires the heat input balancing radiative cooling to be divided between the isentropic and isothermal zones so that each zone is individually balanced. The radiative luminosity L of each zone can be calculated by integrating

$\begin{eqnarray}&&\displaystyle \frac{{dL}}{d\mathrm{ln}r}=4\pi \left(\displaystyle \frac{{n}_{i}}{{n}_{e}}\right){\left(\displaystyle \frac{r}{K}\right)}^{3}{({kT})}^{3}{\rm{\Lambda }}(T)\end{eqnarray} \tag{ 65 }$

over radius. This results in a luminosity

$\begin{eqnarray}&&{L}_{0}=\displaystyle \frac{4\pi }{3}\left(\displaystyle \frac{{n}_{i}}{{n}_{e}}\right){\left(\displaystyle \frac{{r}_{0}}{{K}_{0}}\right)}^{3}{({{kT}}_{0})}^{3}{\rm{\Lambda }}({T}_{0})\,{I}_{3+\lambda }({\alpha }_{K})\end{eqnarray} \tag{ 66 }$

for the isentropic zone and a total luminosity

$\begin{eqnarray}&&L={L}_{0}\left[1+\displaystyle \frac{3\mathrm{ln}({K}_{b}/{K}_{0})}{{I}_{3+\lambda }(1)}\right]\end{eqnarray} \tag{ 67 }$

for the particular case of ${\alpha }_{K}=1$. In order to remain in a configuration with a particular value of K₀, the fraction f₀ of the globally balancing heat input that is deposited into the isentropic zone must be equal to ${L}_{0}/L$.

Specifying a function ${f}_{0}({K}_{0})$ to represent how the feedback required for global balance distributes heat as a function of radius places one more constraint on the system and leads to a potentially unique solution for the global configuration. However, not all functions ${f}_{0}({K}_{0})$ permit a stable solution. Global stability requires a heat-distribution response with the property

$\begin{eqnarray}&&\displaystyle \frac{d\mathrm{ln}{f}_{0}}{d\mathrm{ln}{K}_{0}}\lt \displaystyle \frac{3{f}_{0}}{{I}_{3+\lambda }(1)}\end{eqnarray} \tag{ 68 }$

so that upward perturbations of K₀ lead to net cooling of the isentropic zone and downward perturbations lead to net heating. If that happens, then the system's global configuration remains stable, with a value of K₀ satisfying

$\begin{eqnarray}&&{f}_{0}({K}_{0})\approx {\left[1+\displaystyle \frac{3\mathrm{ln}({K}_{b}/{K}_{0})}{{I}_{3+\lambda }(1)}\right]}^{-1}.\end{eqnarray} \tag{ 69 }$

Qualitatively, the requirement for global stability corresponds to a heat distribution mechanism in which heat deposition onto the isentropic zone does not rise very much, and perhaps even declines, as K₀ rises.

In the context of a configuration with ${\alpha }_{K}\approx 1$, this feature is a physically plausible consequence of feedback energy transport with a bipolar outflow, because the column density of the isentropic zone is ${N}_{0}\propto {r}_{0}{K}_{0}^{-3/2}\propto {K}_{0}^{(2-3{\alpha }_{K})/(2{\alpha }_{K})}$. For systems with ${\alpha }_{K}\gt 2/3$, a rise in K₀ therefore causes N₀ to decrease. The isentropic zone then becomes progressively less able to capture and thermalize the kinetic energy of the outflow as K₀ rises, which enables global regulation.