thermal gas (which the authors called ‘dark Alfven regions’) are found within
the wind, and the transport speed of the CRs is frequently found to lag the
Alfven speed. This echos our claim that the CR flux is, for the most part, not
given by the steady-state form (eqn.8, and see fig.23). The authors attribute
the ‘dark Alfven regions’ to field lines perpendicular to the CR gradient, so
that there is no streaming instability. We note that the low plasma $\beta$ of
the gas implies that the CR acoustic instability (Tsung et al., 2022) can
potentially also play a role, particularly further out in the CGM. This can
also cause regions of flat CR pressure where the CR streaming instability does
not develop, and subsequent onset of the bottleneck effect causes a slowing of
the CRs transport speed (Tsung et al., 2022), which is limited by the Alfven
speed at bottlenecks.

##  6 Conclusions

In this study we explored the effect of CR heating on TI, both in the linear
and nonlinear phase, using a gravitationally stratified setup with vertical
magnetic fields and streaming CRs. In the linear phase, we found that in
accordance to linear theory (Kempski & Quataert, 2020), CR heating can cause
gas entropy modes to propagate at some velocity proportional the Alfven
velocity up the $P_{c}$ gradient. The propagation of the modes is a result of
differential CR heating on different parts of the gas perturbation, resulting
in a net phase velocity. We verified with simulations that the modes propagate
at the expected velocity (see §3.2 and fig.2). This propagation is subdued
when we increase CR diffusivity (fig.3), as CRs diffuse out of the
perturbation before they have time to heat it.

Mode propagation, under the action of CR heating, in a streaming dominated
flow could in principle suppress TI. The idea is that if the time it takes for
the modes to propagate across a cooling radius $t_{\mathrm{cross}}$ is less
than the time it takes the modes to grow, which scales as $t_{\mathrm{cool}}$,
the perturbations would not reach nonlinear amplitudes and become cold clouds.
However, since both the crossing time $t_{\rm cross}\sim L/v_{\rm A}$ and the
heating time $t_{\rm heat}\sim L/v_{\rm A}$ are closely related, there is only
a small range of $\theta=t_{\mathrm{cool}}/t_{\mathrm{heat}}$ near unity where
suppression by propagation effects operates, before the gas becomes over-
heated.

Our most interesting results do not relate to the linear thermal instability,
but rather the non-linear outcome of our simulations. There are two important
things to note. Firstly, thermal instability causes substantial mass dropout
and evolution in CGM properties. A cold disk containing most of the mass forms
in the mid-plane, while the density and pressure of extraplanar gas is
significantly reduced. Thus, initial values of $\alpha=P_{c}/P_{g}$,
$\beta=P_{g}/P_{B}$ evolve substantially during the non-linear stage161616For
instance, in our ‘fountain flow’ simulations, $\beta_{H,i}=300$ becomes
$\beta_{\rm H,f}\approx 1$, due to the reduction in gas pressure. See Fig 17
for other examples.. The strong phase transition from cool disk gas to hot
atmosphere has important consequences for CR wind properties, different from
calculations which assume single-phase isothermal winds. Secondly, the dual
role of CRs in pushing ($\nabla P_{c}$) and heating ($v_{A}\cdot\nabla P_{c}$)
the gas implies that global thermal and hydrostatic equilibrium is not
possible without fine-tuning. Indeed, the gas generally loses force and/or
thermal balance at the disk halo interface, leading to the development of
winds. The sharp reduction in gas density at the disk halo interface
$\rho_{\rm h}\sim\rho_{c}(T_{c}/T_{h})$ reduces radiative cooling
($\propto\rho^{2}$) and gravitational forces ($-\rho g$). The resulting loss
of dynamical and thermal equilibrium can cause gas to accelerate outwards and
heat up. CR winds and fountain flows are most efficient if diffusion is
sufficiently strong such that transport is diffusion dominated until the flow
reaches low densities (i.e., in the halo). Otherwise, since $P_{\rm
c}\propto\rho^{2/3}$ for streaming dominated flows (assuming $B\approx$const),
the sharp density gradient means that CRs suffer strong losses at high
densities, when radiative cooling is still efficient. Since $v_{\rm A}$
increases strongly as the atmosphere thins due to mass dropout, in our
simulations flows are typically streaming dominated in the halo, even if they
are diffusion dominated near the disk.

We find two general classes of solutions, depending on whether the momentum or
heat imparted by CRs dominates. CRs do work by direct acceleration at a rate
$v\cdot\nabla P_{c}$, while the CR heating rate is $v_{A}\cdot\nabla P_{c}$.
Thus, momentum driven winds arise when $M_{A}>1$, and thermally driven winds
arise when $M_{A}<1$. This typically means that momentum or energy driven
winds arise for high and low $\beta$ atmospheres respectively.

  * •

Momentum-Driven Winds; Fountain Flows ($M_{A}>1$; high $\beta$; typically
$\alpha=P_{c}/P_{g}>1$). Cool disk gas is accelerated directly via CR forces –
indeed, winds with almost unchanged characteristics are launched if CR heating
is turned off. Since the flow is super-Alfvenic, magnetic fields are easily
warped; as they wrap around rising gas they can trap CRs, with a consequent
strong jump in CR dominance $\alpha=P_{c}/P_{g}$ once $M_{\rm A}>1$. Density
fluctuations result in a multi-phase fountain flow at low galactic heights,
with the cold dense gas being lifted off the disk and falling back, while
hotter gas flows outwards. At larger distances, when CR heating exceeds
radiative cooling, the gas becomes hotter and mostly single phase. The extent
of the fountain region increases with ${\mathcal{M}}_{\rm A}$ and CR dominance
$\alpha_{\rm H}$.

  * •

Energy-Driven Winds ($M_{A}<1$; low $\beta$; typically
$\alpha=P_{c}/P_{g}<1$). Cool disk gas is strongly heated and evaporated at
the disk halo interface, resulting in a hot wind powered by CR heating. In
steady state, the divergence of the enthalpy flux of the hot gas balances CR
heating. The sharp transition to a hot phase leads to a strong drop in gas
density at the disk halo interface, with $\rho_{\rm h}\sim\rho_{\rm c}(T_{\rm
c}/T_{\rm h})$. This low wind gas density means that the mass flux of CR
heated winds is relatively low; they are an inefficient form of feedback
compared to cool, denser momentum-driven winds, even though a large fraction
of the latter circulates in a fountain flow. The velocity of the wind is
$v\sim 0.4v_{\rm esc}(t_{\rm heat}/t_{\rm ff})^{-1}$, i.e. the flow exceeds
the escape velocity once $t_{\rm heat}\sim t_{\rm ff}$. The flow becomes even
hotter and lower density for these fast winds, leading to very low mass
fluxes. The strong magnetic tension in these sub-Alfvenic flows means there is
little warping of field lines, and the single-phase wind flows monotonically
outward.

There are numerous potential avenues for future work. A key issue is geometry.
We have simulated a plane parallel Cartesian setup. This is appropriate close
to the disk; thus, our simulation domain of $\sim 2H$ is relatively small. As
the flow opens up, a spherical geometry becomes appropriate further away. The
flow properties become quite different, as do the density, velocity, B-field
and hence Alfven speed profiles. Thus, our simulations do not address the
asymptotic properties of the flow far out in the CGM; also, Parker-type sonic
points do not develop in plane-parallel flows. Our work is complementary to 1D
models which use spherical geometry (Ipavich, 1975; Mao & Ostriker, 2018b;
Quataert et al., 2021a, b; Modak et al., 2023), but make other idealizations
which we relax. Other possible extensions include: (i) using more physically
motivated diffusion coefficients which depend on plasma conditions – e.g.,
from quasi-linear self-confinement theory (Wiener et al., 2013a), or using a
model for field-line wandering (Sampson et al., 2022); (ii) incorporating
other sources of thermal and momentum driving besides CRs, and understanding
their mutual interaction; (iii) considering more complex B-field geometry
(e.g., tangled fields due to turbulence). One interesting avenue for future
work include formulating an ‘effective’ 1D model which takes the effect of
multi-phase structure and CR bottlenecks into account, to match our time-
averaged multi-dimensional flows (similar to effective 1D models for turbulent
mixing layers; Tan et al. 2021; Tan & Oh 2021; Chen et al. 2022). Another
would be to make predictions for the nature of CR outflows for different
galaxies lying on the SFR-$M_{*}$ relation (which, in our language, correspond
to different values of $F_{\rm c0}$ and gravity $g$ respectively, leading to
different values of $t_{\rm heat}/t_{\rm ff}$). It would also be interesting
to make observational predictions (e.g., in gamma-ray emission) for CR
dominated fountain flows. Of course, the biggest unknowns are still the
strength of magnetic fields in the CGM, and the nature of CR transport,
particularly the relative importance of streaming and diffusive transport.

## Acknowledgements

We thank Lucia Armillotta, Eliot Quataert, Philipp Kempski, Suoqing Ji, Yan-
Fei Jiang, Eve Ostriker, Ellen Zweibel for helpful discussions. We acknowledge
NSF grant AST-1911198 and NASA grant 19-ATP19-0205 for support. CB was
supported by the National Science Foundation under Grant No. NSF PHY-1748958
and by the Gordon and Betty Moore Foundation through Grant No. GBMF7392.
Additionally, this work made considerable use of the Stampede2 supercomputer
through allocation PHY-210004 from the Advanced Cyberinfrastructure
Coordination Ecosystem: Services & Support (ACCESS) program, which is
supported by National Science Foundation grants #2138259, #2138286, #2138307,
#2137603, and #2138296.

## Data Availability

The data underlying this article will be shared on reasonable request to the
corresponding author. The reader can view videos pertaining to the discussion
in §4 at the following link:
https://www.youtube.com/playlist?list=PLQqhpX30dsYq2cD51L4M2pNQAlm0GSle9.

## References

  * Begelman (1995) Begelman M. C., 1995, in Ferrara A., McKee C. F., Heiles C., Shapiro P. R., eds, Astronomical Society of the Pacific Conference Series Vol. 80, The Physics of the Interstellar Medium and Intergalactic Medium. p. 545
  * Begelman & Zweibel (1994) Begelman M. C., Zweibel E. G., 1994, ApJ, 431, 689
  * Benjamin & Danly (1997) Benjamin R. A., Danly L., 1997, ApJ, 481, 764
  * Booth et al. (2013) Booth C. M., Agertz O., Kravtsov A. V., Gnedin N. Y., 2013, ApJ, 777, L16
  * Breitschwerdt et al. (1991) Breitschwerdt D., McKenzie J. F., Voelk H. J., 1991, A&A, 245, 79
  * Buck et al. (2020) Buck T., Pfrommer C., Pakmor R., Grand R. J. J., Springel V., 2020, MNRAS, 497, 1712
  * Bustard & Zweibel (2021) Bustard C., Zweibel E. G., 2021, ApJ, 913, 106
  * Butsky et al. (2020) Butsky I. S., Fielding D. B., Hayward C. C., Hummels C. B., Quinn T. R., Werk J. K., 2020, ApJ, 903, 77
  * Caprioli & Spitkovsky (2014) Caprioli D., Spitkovsky A., 2014, ApJ, 783, 91
  * Chan et al. (2019a) Chan T. K., Kereš D., Hopkins P. F., Quataert E., Su K. Y., Hayward C. C., Faucher-Giguère C. A., 2019a, MNRAS, 488, 3716
  * Chan et al. (2019b) Chan T. K., Kereš D., Hopkins P. F., Quataert E., Su K. Y., Hayward C. C., Faucher-Giguère C. A., 2019b, MNRAS, 488, 3716
  * Chan et al. (2022) Chan T. K., Kereš D., Gurvich A. B., Hopkins P. F., Trapp C., Ji S., Faucher-Giguère C.-A., 2022, MNRAS, 517, 597
  * Chen et al. (2022) Chen Z., Fielding D. B., Bryan G. L., 2022, arXiv e-prints, p. arXiv:2211.01395
  * Chisholm et al. (2017) Chisholm J., Tremonti C. A., Leitherer C., Chen Y., 2017, MNRAS, 469, 4831
  * Choudhury et al. (2019) Choudhury P. P., Sharma P., Quataert E., 2019, arXiv e-prints, p. arXiv:1901.02903
  * Donahue & Voit (2022) Donahue M., Voit G. M., 2022, Phys. Rep., 973, 1
  * Esmerian et al. (2021) Esmerian C. J., Kravtsov A. V., Hafen Z., Faucher-Giguère C.-A., Quataert E., Stern J., Kereš D., Wetzel A., 2021, MNRAS, 505, 1841
  * Field (1965) Field G. B., 1965, ApJ, 142, 531
  * Fielding & Bryan (2022) Fielding D. B., Bryan G. L., 2022, ApJ, 924, 82
  * Girichidis et al. (2016) Girichidis P., et al., 2016, ApJ, 816, L19
  * Guo & Oh (2008) Guo F., Oh S. P., 2008, MNRAS, 384, 251
  * Guo et al. (2008) Guo F., Oh S. P., Ruszkowski M., 2008, ApJ, 688, 859
  * Hanasz et al. (2013) Hanasz M., Lesch H., Naab T., Gawryszczak A., Kowalik K., Wóltański D., 2013, ApJ, 777, L38
  * Hopkins et al. (2021a) Hopkins P. F., Squire J., Butsky I. S., Ji S., 2021a, arXiv e-prints, p. arXiv:2112.02153
  * Hopkins et al. (2021b) Hopkins P. F., Chan T. K., Ji S., Hummels C. B., Kereš D., Quataert E., Faucher-Giguère C.-A., 2021b, MNRAS, 501, 3640
  * Hopkins et al. (2021c) Hopkins P. F., Chan T. K., Squire J., Quataert E., Ji S., Kereš D., Faucher-Giguère C.-A., 2021c, MNRAS, 501, 3663
  * Hopkins et al. (2021d) Hopkins P. F., Squire J., Chan T. K., Quataert E., Ji S., Kereš D., Faucher-Giguère C.-A., 2021d, MNRAS, 501, 4184
  * Hopkins et al. (2022) Hopkins P. F., Squire J., Butsky I. S., Ji S., 2022, MNRAS, 517, 5413
  * Huang & Davis (2022) Huang X., Davis S. W., 2022, MNRAS, 511, 5125
  * Huang et al. (2022a) Huang T., Stechmann S. N., Torchinsky J. L., 2022a, Phys. Rev. Fluids, 7, 010502
  * Huang et al. (2022b) Huang X., Jiang Y.-f., Davis S. W., 2022b, ApJ, 931, 140
  * Ipavich (1975) Ipavich F. M., 1975, ApJ, 196, 107
  * Jacob & Pfrommer (2017) Jacob S., Pfrommer C., 2017, MNRAS, 467, 1449
  * Ji et al. (2018) Ji S., Oh S. P., McCourt M., 2018, MNRAS, 476, 852
  * Ji et al. (2020) Ji S., et al., 2020, MNRAS, 496, 4221
  * Ji et al. (2021) Ji S., Kereš D., Chan T. K., Stern J., Hummels C. B., Hopkins P. F., Quataert E., Faucher-Giguère C.-A., 2021, MNRAS, 505, 259
  * Jiang & Oh (2018) Jiang Y.-F., Oh S. P., 2018, ApJ, 854, 5
  * Jubelgas et al. (2008) Jubelgas M., Springel V., Enßlin T., Pfrommer C., 2008, A&A, 481, 33
  * Kempski & Quataert (2020) Kempski P., Quataert E., 2020, MNRAS, 493, 1801
  * Kempski & Quataert (2022) Kempski P., Quataert E., 2022, MNRAS, 514, 657
  * Kulsrud & Pearce (1969) Kulsrud R., Pearce W. P., 1969, ApJ, 156, 445
  * Li & Bryan (2020) Li M., Bryan G. L., 2020, ApJ, 890, L30
  * Mao & Ostriker (2018a) Mao S. A., Ostriker E. C., 2018a, ApJ, 854, 89
  * Mao & Ostriker (2018b) Mao S. A., Ostriker E. C., 2018b, ApJ, 854, 89
  * Marasco et al. (2022) Marasco A., Fraternali F., Lehner N., Howk J. C., 2022, MNRAS, 515, 4176
  * McCourt et al. (2012) McCourt M., Sharma P., Quataert E., Parrish I. J., 2012, MNRAS, 419, 3319
  * McCourt et al. (2018) McCourt M., Oh S. P., O’Leary R., Madigan A.-M., 2018, MNRAS, 473, 5407
  * Mertsch (2020) Mertsch P., 2020, Ap&SS, 365, 135
  * Modak et al. (2023) Modak S., Quataert E., Jiang Y.-F., Thompson T. A., 2023, arXiv e-prints, p. arXiv:2302.03701
  * Nelson et al. (2020) Nelson D., et al., 2020, MNRAS, 498, 2391
  * Pakmor et al. (2016) Pakmor R., Pfrommer C., Simpson C. M., Springel V., 2016, ApJ, 824, L30
  * Quataert et al. (2021a) Quataert E., Thompson T. A., Jiang Y.-F., 2021a, arXiv e-prints, p. arXiv:2102.05696
  * Quataert et al. (2021b) Quataert E., Jiang Y.-F., Thompson T. A., 2021b, arXiv e-prints, p. arXiv:2106.08404
  * Quataert et al. (2022a) Quataert E., Jiang F., Thompson T. A., 2022a, MNRAS, 510, 920
  * Quataert et al. (2022b) Quataert E., Thompson T. A., Jiang Y.-F., 2022b, MNRAS, 510, 1184
  * Ruszkowski et al. (2017a) Ruszkowski M., Yang H. Y. K., Zweibel E., 2017a, ApJ, 834, 208
  * Ruszkowski et al. (2017b) Ruszkowski M., Yang H. Y. K., Reynolds C. S., 2017b, ApJ, 844, 13
  * Salem & Bryan (2014) Salem M., Bryan G. L., 2014, MNRAS, 437, 3312
  * Salpeter (1955) Salpeter E. E., 1955, ApJ, 121, 161
  * Sampson et al. (2022) Sampson M. L., Beattie J. R., Krumholz M. R., Crocker R. M., Federrath C., Seta A., 2022, arXiv e-prints, p. arXiv:2205.08174
  * Sharma et al. (2009) Sharma P., Colella P., Martin D. F., 2009, arXiv e-prints, p. arXiv:0909.5426
  * Sharma et al. (2012) Sharma P., McCourt M., Quataert E., Parrish I. J., 2012, MNRAS, 420, 3174
  * Simpson et al. (2016) Simpson C. M., Pakmor R., Marinacci F., Pfrommer C., Springel V., Glover S. C. O., Clark P. C., Smith R. J., 2016, ApJ, 827, L29
  * Skilling (1971) Skilling J., 1971, ApJ, 170, 265
  * Stone et al. (2020) Stone J. M., Tomida K., White C. J., Felker K. G., 2020, ApJS, 249, 4
  * Tan & Oh (2021) Tan B., Oh S. P., 2021, MNRAS, 508, L37
  * Tan et al. (2021) Tan B., Oh S. P., Gronke M., 2021, MNRAS, 502, 3179
  * Thomas & Pfrommer (2019) Thomas T., Pfrommer C., 2019, MNRAS, 485, 2977
  * Thomas et al. (2022) Thomas T., Pfrommer C., Pakmor R., 2022, arXiv e-prints, p. arXiv:2203.12029
  * Tsung et al. (2021) Tsung T. H. N., Oh S. P., Jiang Y.-F., 2021, MNRAS, 506, 3282
  * Tsung et al. (2022) Tsung T. H. N., Oh S. P., Jiang Y.-F., 2022, MNRAS, 513, 4464
  * Uhlig et al. (2012) Uhlig M., Pfrommer C., Sharma M., Nath B. B., Enßlin T. A., Springel V., 2012, MNRAS, 423, 2374
  * Wang et al. (2020) Wang C., Ruszkowski M., Yang H. Y. K., 2020, MNRAS, 493, 4065
  * Wiener et al. (2013a) Wiener J., Oh S. P., Guo F., 2013a, MNRAS, 434, 2209
  * Wiener et al. (2013b) Wiener J., Zweibel E. G., Oh S. P., 2013b, ApJ, 767, 87
  * Wiener et al. (2017a) Wiener J., Oh S. P., Zweibel E. G., 2017a, MNRAS, 467, 646
  * Wiener et al. (2017b) Wiener J., Oh S. P., Zweibel E. G., 2017b, MNRAS, 467, 646
  * Yan & Lazarian (2004) Yan H., Lazarian A., 2004, ApJ, 614, 757
  * Zweibel (2017) Zweibel E. G., 2017, Physics of Plasmas, 24, 055402

##  Appendix A 1D Linearized Equations in Uniform Medium

Kempski & Quataert (2020) showed that 1D calculations can approximately
capture the behavior of CR-modified thermal modes in the linear regime as
perturbations perpendicular to the magnetic field is usually small. In the 1D,
fully-coupled limit, eqn.1-6 reduces to

$\displaystyle\partialderivative{\rho}{t}+v\partialderivative{\rho}{x}+\rho\partialderivative{v}{x}=0,$
(47)
$\displaystyle\partialderivative{v}{t}+v\partialderivative{v}{x}=-\frac{1}{\rho}\partialderivative{P_{g}}{x}-\frac{1}{\rho}\partialderivative{P_{c}}{x}-g,$
(48)
$\displaystyle\partialderivative{P_{g}}{t}+v\partialderivative{P_{g}}{x}+\gamma_{g}P_{g}\partialderivative{v}{x}=-\quantity(\gamma_{g}-1)v_{s}\partialderivative{P_{c}}{x}+\quantity(\gamma_{g}-1)\mathcal{L},$
(49)
$\displaystyle\partialderivative{P_{c}}{t}\quantity(\gamma_{c}-1)\partialderivative{F_{c}}{x}=\quantity(\gamma_{c}-1)\quantity(v+v_{s})\partialderivative{P_{c}}{x}+\quantity(\gamma_{c}-1)\mathcal{Q},$
(50) $\displaystyle
F_{c}=\frac{\gamma_{c}}{\gamma_{c}-1}P_{c}\quantity(v+v_{s})-\frac{\kappa_{\parallel}}{\gamma_{c}-1}\partialderivative{P_{c}}{x},$
(51)

where we have used the fact that the magnetic field $\mathbf{B}$ is constant
in 1D. The streaming velocity
$v_{s}=-v_{A}\mathrm{sgn}\quantity(\partialderivative*{P_{c}}{x})$. Ignoring
contributions from CR sources $\mathcal{Q}$ and assuming the background is
nearly uniform such that we can ignore derivatives of the background but CRs
remain coupled to gas, the linearized equations are

$\displaystyle-\frac{\omega}{\omega_{s}}\frac{\delta\rho}{\rho}+\frac{\delta
v}{c_{s}}=0,$ (52) $\displaystyle-\frac{\omega}{\omega_{s}}\frac{\delta
v}{c_{s}}+\frac{1}{\gamma_{g}}\frac{\delta
P_{g}}{P_{g}}+\frac{1}{\gamma_{g}}\frac{\delta P_{c}}{P_{g}}=0,$ (53)
$\displaystyle\quantity[-i\quantity(\gamma_{g}-1)\frac{\omega_{c}}{\omega_{s}}\Lambda_{T}-\frac{\omega}{\omega_{s}}]\frac{\delta
P_{g}}{P_{g}}-i\quantity(\gamma_{g}-1)\frac{\omega_{c}}{\omega_{s}}\quantity(2-\Lambda_{T})\frac{\delta\rho}{\rho}$
$\displaystyle\quad+\gamma_{g}\frac{\delta
v}{c_{s}}+\quantity(\gamma_{g}-1)\frac{\omega_{A}}{\omega_{s}}\frac{\delta
P_{c}}{P_{g}}=0,$ (54) $\displaystyle\gamma_{c}\alpha\frac{\delta
v}{c_{s}}-\frac{\gamma_{c}\alpha\omega_{A}}{2\omega_{s}}\frac{\delta\rho}{\rho}+\quantity(\frac{\omega_{A}}{\omega_{s}}-i\frac{\omega_{d}}{\omega_{s}}-\frac{\omega}{\omega_{s}})\frac{\delta
P_{c}}{P_{g}}=0,$ (55)

where $c_{s}=\sqrt{\gamma_{g}P_{g}/\rho}$ is the adiabatic sound speed,
$\omega_{s}=kc_{s}$, $\omega_{A}=kv_{A}$,
$\omega_{d}=k^{2}\kappa_{\parallel}$, $\omega_{c}=\rho^{2}\Lambda/P_{g}$,
$\alpha=P_{c}/P_{g}$ and we have assumed the background flow is static.
Substituting $\delta v$ in eqn.52 into eqn.53-55, the equations simplify to

$\displaystyle\frac{\delta P_{g}}{P_{g}}+\frac{\delta
P_{c}}{P_{c}}=\gamma_{g}\quantity(\frac{\omega}{\omega_{s}})^{2}\frac{\delta\rho}{\rho},$
(56) $\displaystyle
i\omega\quantity(\gamma_{g}\frac{\delta\rho}{\rho}-\frac{\delta
P_{g}}{P_{g}})=-\quantity(\gamma_{g}-1)\omega_{c}\quantity[\quantity(2-\Lambda_{T})\frac{\delta\rho}{\rho}+\Lambda_{T}\frac{\delta
P_{g}}{P_{g}}]$
$\displaystyle\quad-i\quantity(\gamma_{g}-1)\omega_{A}\frac{\delta
P_{c}}{P_{g}},$ (57) $\displaystyle\frac{\delta
P_{c}}{P_{g}}\quantity(\omega-\omega_{A}+i\omega_{d})=\gamma_{c}\alpha\quantity(\omega-\frac{\omega_{A}}{2})\frac{\delta\rho}{\rho}.$
(58)

For thermally unstable modes, we expect $\omega\sim\omega_{c}$. If
$\omega_{c}\gg\omega_{s}$, i.e. the cooling time is shorter than the sound
crossing time, $\delta\rho/\rho\ll\delta P_{g}/P_{g}$ from eqn.56, i.e. the
mode is isochoric. If $\omega_{c}\ll\omega_{s}$, the otherwise is true and the
mode is pressure balanced.

The ratio of perturbed CR heating to cooling is

$\frac{\delta\quantity(\mathrm{CRHeating})}{\delta\quantity(\mathrm{Cooling})}=\frac{i\omega_{A}\quantity(\delta
P_{c}/P_{g})}{\omega_{c}\quantity[\quantity(2-\Lambda_{T})\quantity(\delta\rho/\rho)+\Lambda_{T}\quantity(\delta
P_{g}/P_{g})]}.$ (59)

Eqn.59 reduces to eqn.28 in the isochoric and isobaric limits. The ratio
$(\delta P_{c}/P_{g})/(\delta\rho/\rho)$ is given directly by eqn.58

$\frac{\quantity(\delta
P_{c}/P_{g})}{\quantity(\delta\rho/\rho)}=\gamma_{c}\alpha\frac{\omega-\omega_{A}/2}{\omega-\omega_{A}+i\omega_{d}}.$
(60)

For $\omega_{c}\ll\omega_{A}$ (small scale modes) and
$\omega_{d}\gg\omega_{A}$, the RHS is purely imaginary, i.e. $\delta P_{c}$ is
$\pi/2$ out of phase with $\delta\rho$. This phase shift allows CR heating to
counteract cooling.

##  Appendix B Hydrostatic Boundary Conditions for Eulerian Grid Codes

We employ hydrostatic boundary conditions in the $x$-direction, which requires
eqn.23 to be satisfied at the boundaries. In a grid code, cell-center values
are indicated with subscripts $i,j,k$, all of them integers, representing
cells in the $x,y,z$-directions respectively. Below we will shorten the
notation to $i$ to reduce clutter, the relations derived in the following are
implied for all $j,k$. We shall use $is$ and $ie$ to denote the first and last
active zones in the $x$-directions. The cell-center ghost zones are expressed
by $is-n$ and $ie+n$, where $n=1,\ldots,n_{g}$, $n_{g}$ is the number of ghost
zones (typically 2 for piecewise linear method (PLM)). We want eqn.23 to hold
at the boundary cell face, for example at the outer-$x$ boundary

$\evaluated{\derivative{P_{g}}{x}}_{ie+n-1/2}+\evaluated{\derivative{P_{c}}{x}}_{ie+n-1/2}=-\evaluated{\rho}_{ie+n-1/2}\evaluated{g}_{ie+n-1/2},$
(61)

where the fractional index indicates cell faces. The cell-faced values are
approximated linearly as

$\frac{P_{g,ie+n}-P_{g,ie+n-1}}{\Delta
x}+\frac{P_{c,ie+n}-P_{c,ie+n-1}}{\Delta x}\\\
=-\frac{1}{2}\quantity(\rho_{ie+n}+\rho_{ie+n-1})g_{ie+n-1/2},$ (62)

we do not need to approximate $g$ because it is a given function. Using the
ideal gas law $P_{g}=\rho T$ and since our initial profiles are isothermal and
$P_{c}\propto\rho^{\gamma_{c}/2}$ for streaming dominated flows,

$r-1+\alpha\quantity(r^{\gamma_{c}/2}-1)=-\xi\quantity(r+1),$ (63)

where $\alpha,\xi,r$ are defined by

$\alpha=\frac{P_{c,ie+n-1}}{P_{g,ie+n-1}},\quad\xi=\frac{\rho_{ie+n-1}g_{ie+n-1/2}\Delta
x}{2P_{g,ie+n-1}},\quad r=\frac{\rho_{ie+n}}{\rho_{ie+n-1}}.$ (64)

Rearranging,

$\alpha r^{\gamma_{c}/2}+\quantity(1+\xi)r-\quantity(1+\alpha-\xi)=0.$ (65)

Solving for $r$ (e.g. using a non-linear root finder) gives us the value for
$\rho_{ie+n}$ in terms of quantities in the $ie+n-1$ cell-centers, which we
can use to determine $P_{g,ie+n}$ and $P_{c,ie+n}$. $F_{c,ie+n}$ can be
obtained from the CR equation by imposing time steadiness, i.e.

$\evaluated{\derivative{F_{c}}{x}}_{\mathrm{missing}}{bond}=\evaluated{v_{A}}_{\mathrm{missing}}{bond}\evaluated{\derivative{P_{c}}{x}}_{\mathrm{missing}}{bond}$
(66)

Performing a linear approximation,

$F_{c,ie+n}=F_{c,ie+n-1}\\\
+\frac{1}{2}\quantity(v_{A,ie+n}+v_{A,ie+n-1})\quantity(P_{c,ie+n}-P_{c,ie+n-1}).$
(67)

We copy the velocity of the last active zone to the ghost zones if the flow is
outbound and set them to zero otherwise.

The same can be performed for the inner-$x$ boundary, an equation similar to
65 arises,

$\displaystyle\alpha
r^{\gamma_{c}/2}+\quantity(1-\xi)r-\quantity(1+\alpha+\xi)=0,$ (68)

where now $\alpha,\xi,r$ are defined by

$\alpha=\frac{P_{c,is-n+1}}{P_{g,is-n+1}},\quad\xi=\frac{\rho_{is-n+1}g_{is-n+1/2}\Delta
x}{2P_{g,is-n+1}},\quad r=\frac{\rho_{is-n}}{\rho_{is-n+1}}.$ (69)

$F_{c,is-n}$ is given by

$F_{c,is-n}=F_{c,is-n+1}\\\
+\frac{1}{2}\quantity(v_{A,is-n}+v_{A,is-n+1})\quantity(P_{c,is-n}-P_{c,is-n+1}).$
(70)

The magnetic field is set to the value of the nearest $x$-layer. Note that we
assume an isothermal background. If this assumption is relaxed, the energy
eqn.3 will have to be invoked.

Figure 26: Comparison of the fast wind case using different boundary
conditions: hydrostatic (blue) and outflowing (red), as described in the text,
showing no significant difference is found.

In §4 we see that the nonlinear evolution of thermal instability with CR
heating can lead to winds, bringing the system out of hydrostatic equilibrium.
Using the fast wind case as an example, in fig.26 we show that there is no
significant difference when one uses an outflow type condition along the
$x$-boundaries. We copy the density and gas pressure of the last active zone
to the ghost zones and adopted a diode condition for the $x$-velocity (i.e.
the $x$-velocity of the last active zone is copied to the ghost zone if the
gas is outflowing, and zero otherwise). The $y$-velocity at the last active
zone is copied to the ghost zones. As for the CR pressure and $x$-flux, we set
the ghost zone values to be 0.99 times the value at the previous cell. We note
that simply copying the CR pressure and $x$-flux at the last active zone to
the ghost zones would lead to unphysical confinement of CRs as there would be
no CR gradient to transport the CRs out. The CR $y$-flux at the last active
zone is copied to the ghost zones.

##  Appendix C Simulations fixing the base CR flux instead of pressure

As discussed in §4, the nonlinear outcome of TI depends heavily on the CR
heating time $t_{\mathrm{heat}}$ at the halo, which depends on the halo gas,
CR pressure and Alfven speed. In our simulations, we control these quantities
through specifying $\alpha_{0}$, which sets the CR pressure at the base,
$\beta_{0}$, which sets the magnetic field, and $\eta_{H}$, which sets the CR
diffusion coefficient and therefore how much of the base CRs will leak into
the halo. Moreover, these quantities evolve as mass dropout proceeds in the
halo, and the gas pressure $P_{g}$ and Alfven speed evolve. Nonetheless, it is
important to understand the sensitivity of our results to boundary conditions
at the disk. Instead of specifying $\alpha_{0}$, which sets $P_{c}$ at the
base, we could alternatively specify $F_{c0}$, the base CR flux. The
combination of $F_{c0},\beta_{0},\eta_{H}$ will self-consistently determine
what the base $P_{c}$ will be, and is therefore just a different way of
specifying $\alpha_{0},\beta_{0},\eta_{H}$. The advantage of setting
$\alpha_{0},\beta_{0},\eta_{H}$ instead of $F_{c0},\beta_{0},\eta_{H}$ is that
the former involves only dimensionless parameters whereas the latter requires
physical units171717One could specify a parameter like $P_{\rm c,0}/(P_{\rm
g,0}v_{\rm esc})$, but this is similar to our definition of $\alpha_{0}$.. As
we shall demonstrate, our conclusions remain unchanged whether you choose to
fix $F_{c0}$ or $\alpha_{0}$. For example, increasing $F_{c0}$ for a fixed CR
transport model will result in greater presence of CRs in the halo, and is
equivalent to increasing $\alpha_{0}$.

What value of $F_{c0}$ should we set? Let’s assume, in galaxies, all of the
CRs are generated from supernovae. Each supernova releases about $10^{51}\
\mathrm{ergs}$ of energy, for which around $10\%$ goes into accelerating CRs
(Caprioli & Spitkovsky, 2014). Supernovae occur around once per 100 years, so
if we assume all of the CRs produced eventually make it out of the disk, the
rate of CRs released into the halo, would on averaged be $\dot{E}_{c}\sim
10^{51}\ \mathrm{ergs}*0.1/(100\ \mathrm{yr})\sim 3\times 10^{40}\
\mathrm{ergs}\ \mathrm{s}^{-1}$. Assuming CRs escape mostly perpendicular to
the disk, we can relate the CR flux $F_{c}$ with $\dot{E}_{c}$ through
$F_{c}A\approx\dot{E}_{c}$, where $A$ is the galactic disk face area. For a
disk with radius $10\ \mathrm{kpc}$, the CR flux would then be $F_{c}\sim
10^{-5}\ \mathrm{ergs}\ \mathrm{s}^{-1}\ \mathrm{cm}^{-2}$. Let’s convert this
to code units. In our simulations we set $g_{0},T_{0},\rho_{0}$ all to unity
(see §2.2.1). If these variables scale, in real units as $g_{0}=10^{-8}\
\mathrm{cm}\ \mathrm{s}^{-2}$, $T_{0}=10^{6}\ \mathrm{K}$ and
$\rho_{0}=10^{-26}\ \mathrm{g}\ \mathrm{cm}^{-3}$, as typically in galactic
environments, then the scale-height $H=k_{B}T_{0}/m_{u}g_{0}=2.7\
\mathrm{kpc}$, pressure $P_{0}=\rho_{0}k_{B}T_{0}/m_{u}=8.3\times 10^{-13}\
\mathrm{erg}\ \mathrm{cm}^{-3}$, velocity $v_{0}=(k_{B}T_{0}/m_{u})=1=91\
\mathrm{km}\ \mathrm{s}^{-1}$, and the CR flux $F_{c0}=P_{0}v_{0}=1=7.5\times
10^{-6}\ \mathrm{ergs}\ \mathrm{cm}^{-2}\ \mathrm{s}^{-1}$. In code units, the
CR flux $F_{c}\sim 10^{-5}\ \mathrm{ergs}\ \mathrm{s}^{-1}\ \mathrm{cm}^{-2}$
would then be $\sim 1.3$. In reality galaxies could, depending on e.g. the
star formation rate, size and structure, be supplying CRs at different rates,
thus we also explore different values of $F_{c0}$.

Figure 27: Density slice plots at $t=20t_{\mathrm{cool}}$ for a variety of
test cases fixing the base CR flux $F_{c0}$. $F_{c0}$ is given in code-units,
but conversion to real units for scenario specific flows can be found in
Appendix C. The nonlinear outcomes of TI are explored for varying
$F_{c0},\beta_{0},\eta_{H}$. Specifically, the first row displays flows that
are low in $\beta$ and streaming dominated, the second row for low $\beta$ and
higher diffusivity. The third and forth row are replica of the first and
second row at higher $\beta$.

In fig.27 we display the nonlinear outcome of TI for various combinations of
$F_{c0},\beta_{0},\eta_{H}$. Once again we observe the three outcomes
discussed in §4: slow wind, fast wind and fountain flows. The fast wind is
again marked by a single phase, rarefied halo (e.g. middle and right panels of
the second row) while fountain flows are marked by filamentary cold flows
(e.g. right panel of the third row and the middle and right panels of the
bottom row). By varying $F_{c0}$, we can observe the transition into different
outcomes clearly. For example, from a slow wind to a fast wind in the second
row and into a fountain flow in the third and bottom row. All these are the
result of greater supply of CRs to the halo, which increases both CR pressure
support and heating.

Figure 28: $v_{x}$ against $t_{\mathrm{heat}}/t_{\mathrm{ff}}$ at a
scaleheight for the cases shown in fig.27, which the base CR flux is fixed
instead of CR pressure.

In fig.28 we again plot the outflow velocity $v_{x}$ against
$t_{\mathrm{heat}}/t_{\mathrm{ff}}$ (taken at a scale height), recovering the
same trend as in fig.16 that the flow transitions to a fast wind when
$t_{\mathrm{heat}}\ll t_{\mathrm{ff}}$. In short, there is no fundamental
difference whether one fixes $F_{c0}$ or $\alpha_{0}$ at the base. All that
matters is $t_{\mathrm{heat}}$ in the halo.

##  Appendix D Resolution and 3D

Figure 29: Same as fig.16 including cases 2D standard resolution (blue) with
2D higher resolution (red) and lower-resolution 3D (green markings). The
markers indicate ‘slow wind’ (crosses) ‘fast wind’ (diamond) and ‘fountain
flow’ (triangular) profile parameters respectively. Test cases with flow
parameters other than the three mentioned are indicated by grey circles and
are shown for reference. Figure 30: Time averaged projection plots of the
density, velocity and temperature for the ‘slow wind’, ‘fast wind’ and
‘fountain flow’ cases. In each case, the fiducial profile (blue) is compared
against the higher resolution profile (red) and the 3D profile (green).

We rerun the ‘slow wind’, ‘fast wind’ and ‘fountain flow’ cases in §4 with
higher resolution and in 3D to check that our results hold. For increased
resolution, we resolve the simulation domain ($-2H<x<2H$) by $2048\times 512$
grids (doubling the $x$-resolution) whereas for 3D simulations, the grid
resolution is reduced to $256\times 128\times 128$ (again higher resolution
along the $x$-axis) to save computational time. Further details regarding the
setup are listed in table 1. As shown in fig.29, the fluid properties of the
higher resolution and 3D runs are all in line with the trends given by
simulations with fiducial resolution. Our conclusions remain unchanged. While
there is some scatter as resolution and dimensionality change, the scatter
lies along the trends we have already found. In fig.30 we compare the time
averaged projection plots of the density, velocity and temperature for all
three solution outcomes. The fiducial and the high resolution 2D profiles are
very similar. The 3D profiles also give very similar outcomes. The time-
averaged profiles do deviate somewhat more: for instance, up to a factor of
$\sim 2$ in asymptotic temperatures or densities, though at least part of this
is due to the much lower resolution in our 3D sims.