A Review of the Theory of Galactic Winds Driven by Stellar Feedback

Galactic winds from star-forming galaxies are crucial to the process of galaxy formation and evolution, regulating star formation, shaping the stellar mass function and the mass-metallicity relation, and enriching the intergalactic medium with metals. Galactic winds associated with stellar feedback may be driven by overlapping supernova explosions, radiation pressure of starlight on dust grains, and cosmic rays. Galactic winds are multiphase, the growing observations of emission and absorption of cold molecular, cool atomic, ionized warm and hot outflowing gas in a large number of galaxies have not been completely understood. In this review article, I summarize the possible mechanisms associated with stars to launch galactic winds, and review the multidimensional hydrodynamic, radiation hydrodynamic and magnetohydrodynamic simulations of winds based on various algorithms. I also briefly discuss the theoretical challenges and possible future research directions.


Introduction
The basic scenario of galaxy formation and evolution in the cosmological context can be described in the framework of the standard Λ cold dark matter (ΛCDM) cosmology. However, many observed properties of galaxies remain difficult to explain. To obtain a realistic picture of cosmological galaxy formation and evolution which matches observations, astrophysical processes such as gas cooling, accretion, galaxy mergers, and feedback should be considered. Among them feedback is one of the most difficult ingredients to be understood. It is widely believed that feedback provides energy, mass, and metal to the interstellar medium (ISM) and the circumgalactic medium (CGM), and drives galactic-scale winds.
Feedback from both stars and active galactic nucleus (AGN) activities has been proposed to explain galactic winds. The discovery of galactic winds can be traced back more than a half century ago, when a large-scale outflow in a starburst galaxy M82 was observed [1,2]. At almost the same time quasars were discovered [3]. Although the energy source of quasars was heavily debated in the beginning, the theory that quasars and other AGNs are powered by gas accretion onto supermassive black holes (SMBHs) was finally accepted by the end of 1960s [4][5][6][7][8][9]. Since then, galactic winds have been observed ubiquitously in both the local and high-redshift Universe [10][11][12]. The theory for galactic winds developed over years is of outflows powered by the energy and momentum of multiple supernovae (SNe) from massive stars [13][14][15][16][17][18][19]. By the end of the last century, AGN feedback and AGN-driven outflows have also gained momentum in the literature [20][21][22][23][24].
AGN feedback is the astrophysical process to link the momentum and energy released by the central SMBHs to the surrounding medium, and drives outflows (see reviews [25][26][27][28] ). The energy

Impact of Galactic Winds Driven by Stellar Feedback
Galactic winds driven by stellar feedback are important on various scales. On cosmological scales, the most significant effect of galactic winds is to shape the cosmic galaxy luminosity function, flattering its low-mass end slope compared to that of the halo mass function. The observed galaxy luminosity function shows a flat low-mass-end slope and sharp exponential cut-off at the high-mass-end, which are different from the shape of the dark matter halo mass function predicted by the ΛCDM model [59][60][61]. The high-mass-end can be explained by the AGN feedback [33,62,63], while the flat low-mass-end of the galaxy luminosity is believed to be caused by stellar feedback, which drives gas out of galaxies, and suppress the formation of low-mass galaxies [13,59,[64][65][66][67].
Also, galactic winds play a crucial role in determining the chemical evolution of galaxies and mass-metallicity relation. Galactic winds from stellar feedback are expected to be metal-rich, and massive galaxies with deeper gravitational potential are expected to retain more of the galactic wind ejecta than dwarf galaxies, therefore metal loss from galaxies is strongly anti-correlated with the mass of galaxies [60,[67][68][69][70][71][72][73]. Galactic winds are also suggested to be the primary sources to pollute the circumgalactic medium (CGM) and the intergalactic medium (IGM) with metals [74][75][76][77].
Galactic winds are probably the most extreme manifestation of the feedback between star formation in a galaxy and its ISM [78][79][80][81]. On galactic scales, winds can heat gas in the ISM and prevent gas from cooling and forming stars in both dwarf and massive galaxies [19,82]. The unbound outflowing gas is injected into the CGM and IGM. On the other hand, if the velocity of a wind is below the galactic escape velocity, the wind eventually falls back onto the galaxy to form galactic fountain or later gas accretion inflows, which helps to establish the so-called "gas recycling" process in galaxies [83][84][85][86]. Since galactic winds are from the central region of galaxies with low angular momentum, the recycling of gas which falls back from the galactic inner region to outer region re-distributes angular momentum in the galaxies [86][87][88].
X-ray satellites including Chandra and XMM-Newton observatories have been used to observe highly ionized X-ray emitting hot phase of galactic winds, and also to constrain the hard X-ray emission directly from very hot winds driven by supernova (SN) explosions [89][90][91][92][93]. Emission lines such as Hα, N II, O II, O III, and absorptions such as Na I D, K I, Mg II have been used to trace morphologies and velocities of warm, cool and cold outflows in galactic winds, including warm ionized [94][95][96][97][98][99][100][101][102], neutral atomic [19,[103][104][105][106][107][108][109][110], molecular [111][112][113][114][115][116][117][118], and dust [119][120][121][122][123][124] outflows. As shown in the references, the multiphase outflowing gas can be accelerated to hundreds of km s −1 , even above 1000 km s −1 . Taking the starburst M82 as an example: Figure 1 shows the Hα traced warm phase of galactic wind in M82, and Figure 2 demonstrates the schematic of the multiphase wind around M82. The Hα emitting clouds reach ∼600 km s −1 [94,125], and the cool and cold outflows have velocities ∼100-300 km s −1 e.g., [117,126]. Also, the galactic wind in another starburst NGC 253 show similar properties [114,118]. A review of the recent observations of multiphase stellar-driven galactic winds by David Rupke can be found in this special Issue. Previous research suggested that the warm phase of the wind is accelerated by the ram pressure of hot phase driven by SN explosions e.g., [19,92], but recent simulations showed that such ram-pressure-acceleration scenario may be problematic [127][128][129][130][131]. From a theoretical point of view, the driving mechanisms of multiphase winds driven by stellar feedback are far from clear, and the comparison between theory and observations is still incomplete.  observational evidence for the red statements. Similar schematic has been made for NGC 253 [113,116]. Below I give a brief theoretical review of galactic winds driven by stellar feedback. In Section 2, I review the current models of very hot galactic winds driven by SN explosions, the interaction between hot winds and the multiphase ISM, and the entrainment and acceleration of clouds in hot winds. In Section 3, I discuss the possibility of driving large-scale dusty winds by radiation pressure produced by the continuum absorption and scattering of starlight on dust grains. In addition to analytic models, I also show the radiation hydrodynamic simulations of the coupling between radiation, dust, and gas, which results depend on the numerical algorithms. In Section 4, I introduce cosmic rays as an alternative mechanism to launch and propagate galactic winds. Although the physics of cosmic-ray transport is still unclear, many theoretical models and numerical simulations have been carried out to explore cosmic-ray feedback in various galaxies. Conclusions and discussion are in Section 5.

Impact of Supernovae and Supernova Remnants on the ISM
The best-developed model of galactic winds invokes the heating of the ISM by overlapping core-collapse SN explosions. Before I discuss the impact of net SN explosions, I first briefly review the dynamics and feedback of an individual supernova remnant (SNR). A SN from a massive star typically ejects 2-10 M material into the ambient ISM and drive a shock into the ambient ISM. The interaction between a SN and its surrounding medium forms a SNR which can be observed by multiwavelength telescopes. The spherical expansion of a SNR in a uniform medium has been systematically studied and is divided into several stages: the free expansion, Sedov-Taylor, pressure-driven snowplow, and momentum-conserving stages (see review by Ostriker & McKee [132]).
During the initial free expansion stage, the SNR expanses with nearly constant velocity until the swept-up mass of the ambient medium is comparable to the ejecta mass. Then the SNR enters the second stage-the adiabatic Sedov-Taylor stage, which was first studied by Sedov [133] and Taylor [134].
In this stage the information of SN ejecta mass is erased, and the properties of the SNR mainly depends on the ambient density and the initial SN energy. The structure of the SNR can be described by a self-similar solution. Let the radius of the forward shock front to be r s , the self-similar structure gives r s ∝ E α n β t η where E, n, t are the initial SN energy, ambient density and time, and α, β, η are the self-similar indexes which can be solved by dimensional analysis. The solutions are [135,136] where E 51 = E/10 51 ergs, t 3 = t/10 3 yr and n 0 = n/cm 3 . Once the radiative cooling in the SNR becomes important, the SNR leaves the Sedov-Taylor stage. According to Kim & Ostriker [136], the cooling time t cool , the radius of the forward shock front r cool , the shock front velocity, the post-shock temperature T cool , and the swept-up medium mass M swept are given by T cool = 5.67 × 10 5 K E 0.13 51 n 0.26 0 , For t > t cool , the SNR goes to the pressure-driven snowplow stage, in which the thermal energy of the SNR interior still cannot be neglected. The expansion of the SNR is due to thermal pressure of the SNR interior region, and the radius and velocity of the SNR shell can be approximately given by [135] r slowplow ≈ r cool (t/t cool ) 2/7 , The SNR leaves the pressure-driven snowplow stage once the thermal energy of the SNR is radiated away, but the SNR still expands until its forward shock speed approaches the effective sound speed in the ambient medium (c s ). The SNR fade-away time and size can be estimated by [137] t fade ≈ 1.83 × 10 6 yr E 0.32 51 n −0.37 where c s,10 = c s /10 km/s. The analytic estimate of SNR evolution in a uniform medium has been tested and confirmed by one-and multidimensional hydrodynamic simulations of spherical SNRs with radiative cooling [136,[138][139][140][141][142][143]. However, since the realistic ISM is highly inhomogeneous [144], it is important to investigate the expansion of SNRs in an inhomogeneous medium [137,142,143,145]. Martizzi et al. [142] initialized the turbulent medium with a lognormal density distribution [146,147] and performed 3D hydrodynamic simulations with adaptive mesh refinement treatment to compare SNR expansions in an inhomogeneous medium to that in a uniform medium. They found that no clear stages can be divided for SNRs in an inhomogeneous medium. At the same mean density, SNRs in an inhomogeneous medium expand faster along some preferred directions due to the existence of low-density channels compared to those in a uniform medium, and the momentum feedback provided by SNRs in an inhomogeneous medium is ∼70% of that in a uniform medium. Zhang & Chevalier [137], on the other hand, did not observed low-density channels in an inhomogeneous medium, but found that the properties of remnants are mainly controlled by the mean density of the ambient medium.
Similar to galactic winds, the ISM is referred to have multiphase structure. Without SN feedback, the ISM medium can be separated into two phases: the cold neutral medium and warm medium as the consequence of thermal instability, while the warm phase can be further divided into warm neutral and warm ionized medium [148][149][150][151][152]. The two-phase ISM model was later extended to include the feedback from SN explosions. Cox & Smith [153] first proposed that SNRs may create large-scale "tunnel" which is filled with hot and tenuous medium. McKee & Ostriker [154] then proposed the famous three-phase ISM model. The third phase is the "hot ionized medium", which had been shock heated by overlapping SNRs to ∼10 6 K and occupies most of the volume of the ISM. According to the three-phase model, the warm and cold medium is embedded in the hot medium as clouds, which are supposed to be in thermal pressure equilibrium with the hot medium. However, observations have shown that the hot medium may be highly overpressured compared to other phases [155], thus the interstellar magnetic fields may play an important role in equilibrating the pressure balance and suppressing the thermal evaporation of warm and cold clouds. Another problem of the three-phase model is that it cannot explain the substantial amount of warm H I in the ISM [135]. Although many details are still under debate (see discussion by Cox [156]), the three-phase ISM model has been the dominant conceptual framework of the ISM for four decades, and recent numerical simulations have provided better understanding of the multiphase ISM [157][158][159][160]. The detailed discussion of the ISM model is beyond the scope of this review, but keep in mind that the theory of galactic winds is based on the theory of the ISM.

The CC85 model
According to McKee & Ostriker [154], a galactic-scale "superbubble" forms if the SNR net energy cannot be efficiently radiated away by the multiphase ISM, otherwise a "galactic fountain" is produced to interact with the galactic halo. Although the physical processes of galactic wind launching and propagating by SN explosions are complex, a simple, clear analytic estimate can be applied to galactic winds. The classical SN-driven galactic wind model was first developed by Chevalier & Clegg [161] (hereafter CC85 model).
Assuming a spherical symmetry with negligible radiative cooling, the total mass and energy injection rate into the wind areṀ hot andĖ hot respectively, the steady-state hydrodynamic equations for a hot wind driven by overlapping SN explosions are 1 where r, u, ρ, P are the wind velocity, radial radius, density, and pressure, respectively. For r < R the averaged injected efficiencies per unit volume q and Q are defined as q =Ṁ hot /V and Q =Ė hot /V, where V = 4πR 3 /3 is the volume of the starburst/wind launching region with R being the radius of this region. For r > R, we set q = Q = 0. The solutions of Equations (13)- (15) can be translated to the equations of the Mach number for the wind M = u/ γP/ρ, which satisfies for r < R and for r ≥ R. The solution of M smoothly transfers from subsonic flow M < 1 for r < R to supersonic flow M > 1 for r > R. The critical point M = 1 is at r = R. The transonic solution is similar to the Parker solar wind model, in which the velocity of the wind increases from subsonic to supersonic at a critical distance [162]. However, in contrast to the Parker's model which is in an isothermal process, the CC85 model applies to adiabatic flows, while the wind structure changes for radiative flows (see Section 2.3). The wind solution is a function of radius r. We take r * = r/R as the dimensionless radius, and dimensionless velocity u * , density ρ * and pressure P * as the solutions of these dimensionless variables are shown in Figure 3. For r R (r * 1), we have u * ∝ r * , and for r R (r * 1), the asymptotic velocity has u * → √ 2, ρ * → 1/(4 √ 2πr 2 * ) ∝ r −2 * , and P * ∝ r −2γ * .  Next, we introduce two dimensionless parameters for the CC85 model, the thermalization efficiency α and the mass-loading rate β asĖ where SFR is the total star formation rate in a galaxy, the net energy rate provided by SNeĖ SN can be estimated by SFR thatĖ where e = 10 51 e 51 ergs being the energy injected by an individual SN, and ν = (100M ) −1 ν 100 is the number of SNe per unit mass of star formation. Typically, for 100 M of star formation one SN explodes, thus ν 100 1. More specifically, ν 100 depends on the stellar initial mass function (IMF) and we have ν 100 = 1.18 for Salpeter IMP [163]. We define α = α e 51 ν 100 (24) as the new thermalization efficiency, thus we have Therefore, the asymptotic velocity of the wind is V ∞ hot = 10 3 km s −1 α 1/2 β −1/2 , the density and temperature of the wind at the center of the galaxy are n hot = 3 cm −3 α −1/2 β 3/2 R −2 200 pc SFR 1 and T hot = 1.5 × 10 7 K(α/β), where SFR 1 =SFR/10 M yr −1 . Please note that for α ∼ β, the wind emits hard X-rays, which are hotter than the superbubble generated by the hot ionized medium according to the three-phase ISM model.
In addition to the CC85 model on galaxy scales, Cantó, Raga & Rodríguez [164] developed an analytic model for winds from star clusters, where the multiple stellar winds interact with clusters of massive stars. Although the sources and scales of the cluster winds are different from the SN-driven galactic winds, the model is essentially the same as the CC85 model. Numerically, Strickland & Heckman [92] performed a series of hydrodynamic simulations in 2D cylindrical symmetry to model the SN-driven winds from disk galaxies and compared to the results of the spherical CC85 model. They found that the wind from a disk-like starburst region of diameter 750 pc and height 105 pc has equivalent structure to that of a spherical CC85 model. The hot gas properties within a disk-like starburst region can be well predicted using the CC85 model. Also, recent 3D global simulations of SN-driven winds showed that the properties of SN-driven winds in the adiabatic phase are consistent with the CC85 model [165,166]. Therefore, it is widely believed that the CC85 model is a good approximation for adiabatic winds.
Importantly, there are two controlling parameters in the CC85 model: the thermalization efficiency α, and the mass-loading rate β, both of which are difficult to be determined by direct observations. One case study is for M82. Strickland & Heckman [92] used Chandra X-ray observatory to constrain the hot wind from M82. The diffuse hard X-ray directly shows the existence of hot flows with the temperature of ∼3-8 × 10 7 K, the thermalization efficiency and mass-loading rate can be constrained as α ∼ 0.3-1 and β ∼ 0.2-0.6. More generally, Zhang et al. [167] derived a general constraint on β across a wide range of galaxies from dwarf starbursts to luminous infrared galaxies (LIRGs) and ultraluminous infrared galaxies (ULIRGs) using the observed linear relation between the X-ray luminosity and SFR [168][169][170][171][172][173][174]. They found that for galaxies with SFR 10 M yr −1 the mass-loading rate is β 1. This constraint means that the efficiency of converting the ambient ISM into the material of the hot wind fluid by overlapping SNe is low according to the CC85 model, which alone cannot explain the high mass-loading β ∼ 1-10 required by the integrated constraints on stellar feedback models in a cosmological context [65,69,[175][176][177]. Bustard et al. [178] modified the CC85 model including non-uniform mass and energy distribution in the starburst region, and radiative cooling in hot flows. They also used the L X −SFR relation to constrain β and found that β 5 for SFR 10 M yr −1 in general, and β 2 for SFR 10 M yr −1 if the X-ray emission contributed by hot wind is less than 10% of the total X-ray emission. Please note that the above constraint on β is only for mass-loading rate of the very hot phase of galactic winds (hereafter β hot ). The constraint on mass-loading rates of other wind phases can be directly estimated by observations [96,99,100,102,104,105]. For example, Martin [104] estimated the mass-loading rate of neutral atomic outflows from ULIRGs with SFRs of 100-1000 M yr −1 , and found that the mass-loading rate β cool ∼ few, which is one order of magnitude higher than β hot . Rubin et al. [99] found a more conservative estimate of the mass-loading rate of warm phase β warm 0.02-0.6 from star-forming galaxies with SFRs 100 M yr −1 at redshift z ∼ 0.5, while Weiner [96] found that the mass-loading of warm ionized phase β warm ∼ 1 for star-forming galaxies at z ∼ 1. Heckman et al. [100] had a higher estimate for the warm ionized phase β warm ∼ 1-4 from starburst galaxies. Chisholm et al. [102] found the warm phase of winds has scaling relations between mass-loading rates and the stellar mass of their host galaxies β warm ∝ M −0.4 * . This scaling relation is similar to recent numerical simulations, but the scalings obtained by observations are a factor of five smaller than some simulations [81]. In general, a large fraction of the mass of galactic winds are in other phases rather than the very hot phase. Then a question arises: what are the launching and acceleration mechanisms for the multiphase winds?
The mass-loading rate and thermalization efficiency of hot wind have also be constrained by hydrodynamic simulations [165,[179][180][181][182][183][184][185]. Some work found low mass-loading rates β 1 but others found higher. Usually these hydrodynamic simulations combined the very hot and hot components and defined hot winds with temperature a few ×10 5 K. Many simulations used local Cartesian boxes with vertical stratification for the ISM structure and injected SN explosions in the boxes to model the SN-driven winds. Creasey et al. [179]'s simulations explored mass-loading rates as a function of galactic disk surface density and found that β β is for all phases of winds, but not only for hot winds. Martizzi et al. [181] performed local simulations and found low mass-loading rates β 1. However, they showed that the local boxes cannot well predict the properties of hot winds because they do not capture the realistic global geometry and gravitational potential of galaxies. Fielding et al. [165] used a global galactic disk setup with more realistic gravity and geometry than the vertically stratified structure to revisit the problems in Martizzi et al., and found that α 10 −2 and β 1, and α and β decrease with the surface density Σ gas . However, β decreases with spatial resolution while α is well converged with resolution. Li et al. [184] performed simulations with fixed SN rates and more sophisticated treatment of radiative cooling, and found β hot 2.1(Σ gas /1 M pc −2 ) −0.61±0.03 with total mass-loading rate β 7.4(Σ gas /1 M pc −2 ) −0.61±0.03 , which shows a shallower decline than that in Creasey et al. [179]. The hot wind in their simulations is defined with flow temperature T > 3 × 10 5 K, then, for solar neighborhood Li et al.'s model gives β hot 0.8 and β 2-3. Kim & Ostriker [185] gave a lower mass-loading rate with the averaged β ∼ 0.1 for solar neighborhood. They defined hot winds as flows with temperature T > 5 × 10 5 K. Overall, the above numerical results basically support the constraint β hot 1. On the other hand, higher value of β was also obtained. Gatto et al. [182] defined hot winds with temperature T > 3 × 10 5 K, and found that thermal pressure-driven outflows have β 1 at 1 kpc if half of the volume near the galactic disk mid-plane can be heated to the hot phase. Girichidis et al. [180] simulated winds for solar neighborhood condition and found large mass-loading rates β ∼ 5-10. According to Kim & Ostriker [185], the reason for these large differences in mass and energy loading among many authors are mainly due to the different simulation settings in the vertical scale height of SNe.

The CC85 Model with Radiative Cooling
The traditional CC85 model is for adiabatic flows. Assuming the radius of a hot wind is much larger than the radius of wind launching/starburst region r R, the dynamic timescale of a hot wind at radius r is given by where r 10 = r/10 kpc. On the other hand, using an approximated power-law for radiative cooling, the radiative cooling timescale of the wind at r is (see Thompson et al. [186] for details) Remembering that the CC85 model is for a spherical wind, but Equation (29) includes the factor of opening angle of the wind Ω 4π = Ω/4π as a generalized case. Setting t cool = t adv , we can derive the critical radius beyond which radiative cooling is important in the flow: Equation (30) shows that the cooling radius sensitively depends on both α and β. Setting r cool = R, we can also derive a critical β which is given by For example, for ULIRGs with SFR = 100 M yr −1 and R = 1 kpc, we have β crit ∼ 0.28(α/0.5) 0.636 . For M82-like starburst with SFR ≈ 10 M yr −1 and R = 300 pc, β crit ∼ 0.41(α/0.5) 0.636 . If β β crit , the hot flow is radiative outside the starburst region R, and the CC85 model should be modified.
In contrast to the CC85 model, no analytic solutions can be obtained for the wind model with radiative cooling. Semi-analytic and numerical methods have been used to study the impact of radiative cooling on galactic or cluster winds, first by Wang et al. [187,188], then by others [178,186,[189][190][191][192][193]. Compared to the solution of the CC85 model, some significant differences are as follows. The temperature of a wind drops quickly from ∼10 7 K to ∼10 3 -10 4 K due to rapid cooling, therefore the sound speed drops, and the Mach number increases. The sonic point is no longer at the starburst radius R. Since the hot wind transits from the very hot phase to the warm phase, Thompson et al. [186] developed a model that the warm clouds in galactic winds are produced by the cooled hot winds due to thermal instability (see Section 2.4 for more details). The X-ray luminosity is also dimmer than that without radiative cooling, therefore the X-ray observations gave a less stringent constraint on β hot than that from the CC85 model [178].
In addition to radiative cooling, Scannapieco [194] found that Compton cooling due to scattering between free elections and photons from the host galaxy can be more important than radiative cooling by atoms and ions in hot winds. The atomic and ionic radiative cooling can enhance density inhomogeneities in hot winds. Therefore, although the main cooling mechanism is Compton cooling, the hot wind can eventually generate warm clouds with a temperature of ∼10 4 K, similar to the scenario given by Thompson et al. [186].

Acceleration of Warm, Cool and Cold Clouds in Hot Winds
Galactic winds have been observed to be multiphase, with clear evidence for cool, warm, and cold components from multiwavelength observations. A key question is how to accelerate the multiphase gas to the observed velocities with a typical value of a few hundred km/s, or even a thousand km/s. Since a hot wind driven by SNe is highly supersonic, the ram pressure of the wind is much stronger than the thermal pressure of it. A prevailing scenario is that the ISM material is advected into a hot wind driven by SNe and accelerated by the ram pressure of the wind [19,89,92,104,[195][196][197].
Global simulations on galaxy scales have been carried out to model the acceleration of the ISM by a hot wind driven from the starburst region. In particular, the Hα emitting ionized warm and Na I D cool phase of winds were modeled recently. Cooper et al. [196] performed 3D simulations to study the interaction between a hot wind and an inhomogeneous ISM. They found that Hα filaments are generated due to the breakup of clouds in the ISM in starburst region and accelerated by the ram pressure of the wind. In addition to the mass-loaded hot wind with big β, they proposed that soft X-ray also arise from bow shocks upstream of the dense clouds advected into the hot flow, and the interaction between these shocks. Cooper et al.'s simulations showed that the soft X-ray emission region is basically associated with the Hα emitting gas. Compared to the observation, the Hα emitting clumps with a velocity of ∼600 km s −1 in M82 was considered to be accelerated by the hot wind [92]. Moreover, Fujita et al. [197] carried out 2D simulations to study the origin of the Na I D-absorbing gas in ULIRGs. They found that cool gas can be accelerated to a mean velocity of ∼320 km s −1 , and some gas can be accelerated to 750 km s −1 . Therefore, they claimed that their model can explain the kinematics of cool gas seen in the Na I D observations. It seems that global simulations as a powerful tool have demonstrated the acceleration of a multiphase wind and matched observations. However, high-resolution simulations of cloud-hot wind interaction discovered a different story.
Motivated by global simulations which showed the formation and acceleration of filaments and clumps in hot winds, the dynamics of an individual cloud entrained in a hot wind has been studied. A more generalized problem-the interaction between clouds and hot medium or shocks in the context of ISM and hot winds have been studied for decades, both analytically [129,[198][199][200], and numerically [127,128,131,[201][202][203][204][205][206][207][208][209][210][211][212][213][214][215][216][217][218][219]. Although many works suggested that a cloud can survive in hot flows or shocks, recent numerical simulations focus on clouds on pc or sub-pc scales and hot flow parameters of galactic winds including radiative cooling and thermal conduction showed that a cloud is quickly shredded by the hot wind [127,128,131]. To describe the dynamics of clouds, an important timescale needs to be introduced is the cloud crushing time t cc , which is the characteristic timescale for the cloud to be crushed by the hot wind passing through the cloud where R c and ρ c are the radius and density of the cloud. Klein et al. [200] showed that the cloud destruction timescale is a few t cc . In the meantime, the cloud is being accelerated by the ambient hot wind. The timescale for the cloud to be comoving with the hot wind, i.e., the acceleration timescale is For ρ c ρ hot , we always have t acc t cc . If the destruction time of a cloud is comparable to t cc , the cloud can never be fully accelerated by the hot wind. However, numerical simulations showed that the cloud lifetime can be prolonged due to radiative cooling, thermal conduction, and magnetic fields.
Some more recent work suggested that clouds entrained in a hot wind are disrupted by the Kelvin-Helmholtz (KH) instability on a timescale of t KH ≈ κt cc , with κ ∼ 10 when radiative cooling in the cloud is efficient [211,220]. Scannapieco & Brüggen [127] used the cloud-following scheme to study the evolution of cold clouds entrained in a hot wind, and found that the KH instability is strongly damped in supersonic hot flows, while clouds are shredded by shearing instability, which has a timescale ∝ t cc √ 1 + M hot with M hot being the Mach number of the hot wind. In particular, they found that the timescale for half of a cloud to be below 2/3 of the initial cloud density is t half ≈ 4t cc √ 1 + M hot , which is the typical lifetime of the cloud. Schneider & Robertson [130] did similar simulations for both turbulent and spherical clouds but found a longer lifetime for spherical clouds. The difference is due to the different treatment of radiative cooling between two groups. Moreover, thermal conduction may also be important to evaporate the clouds [128,198,221]. Brüggen & Scannapieco [128] found that the electron thermal conduction can extend the lifetime of clouds by compressing clouds into dense filaments. They discussed that since the compression by thermal conduction causes the clouds to present a small cross section to hot flow, the terminal velocity of the clouds is still far below the hot wind velocity. Ignoring thermal conduction, McCourt et al. [131] also observed the fragmentation of clouds interacting with a hot wind but argued that these fragments increase the total surface area and may cause rapid entrainment of cold gas. However, the fragments with a typical column density of ∼10 17 cm −2 is too small to be resolved by their simulations. Whether such small fragments can be accelerated to high velocity is still an open question.
Some studies have shown that magnetic fields can delay the destruction of clouds [202,[214][215][216]218]. In addition to the ram pressure of hot winds, McCourt et al. [215] showed that the magnetic field in a hot wind with M hot = 1.5 enhances the drag force on an entrained cloud by a factor of ∼(1 Banda-Barragán et al. [216,218] developed a model of magnetized clouds entrained in a hot flow with higher Mach number M hot 4-5, and found that spherical clouds may be accelerated to ∼0.1-0.2 of the hot wind velocity for the plasma beta β A 10-100, which corresponds to a magnetic field strength of ∼3-10 µG α 1/4 β 1/4 R −1 200 pc SFR 1 1/2 . Turbulence also helps the cloud expansion and increases the cross section of the cloud, leading the cloud to reach ∼0.3-0.4 of the hot wind velocity. Magnetic fields seem to be the key ingredient to accelerate the ISM clouds in a hot wind, and the realistic magnetic field structure constrained by observations need to be used for magnetohydrodynamic (MHD) simulations (e.g., M82 magnetic fields [222]). Please note that even magnetic fields are important to suppress the shredding and evaporation of the clouds, the spatial extent and acceleration profile should still be tested and compared to observations of resolved systems such as M82 and NGC 253 [117,118]. Another possible origin of warm/cool gas in galactic winds arises from thermal instability in the winds. Returning to the radiative cooling of hot winds (Section 2.3), Thompson et al. [186] proposed that the multiphase gas from the ISM is initially accelerated by the ram pressure of a hot flow, then rapidly shredded, thereby increasing the mass-loading of the hot flow, and then the thermally driven hot flow with radiative cooling produces warm/cool gas again with high velocities by thermal instability. The cooled wind interacts with the CGM and deposit warm/cool gas into the halo of the host galaxy at ∼50-100 kpc. This result may explain the observed warm/cool outflows in the halos of galaxies [223][224][225][226][227][228]. Figure 4 shows the schematic of the evolution of a hot wind with radiative cooling. This picture was studied numerically by Schneider et al. [166], who found multiphase wind generated in their simulations. However, note that this picture assumes that the ISM near the starburst/wind launching region is destroyed completely and only contributes to the mass-loading of the wind. It is important to explore the interaction between the ISM and the hot wind on galaxy scales more carefully in the future. As mentioned in Section 2.3, Scannapieco et al. [194] observed inhomogeneities and fragments formation in hot flows due to Compton and atomic/ionic cooling. This picture may provide another explanation of warm/cool gas in the halos, but also neglects the interaction between the wind and the ISM and CGM.

Radiation Feedback, and Radiation Pressure on Dust
Radiation from individual or multiple stars may also be important in regulating star formation and launching galactic winds.
For massive stars, photoionization and radiation pressure are the dominant feedback mechanisms to affect their surrounding medium before the final SN explosions. Young massive stars emit mostly in  . The schematic of the evolution of a hot wind to a large distance of ∼100 kpc. The hot wind launched from the starburst region (host) cools radiatively above the cooling radius at ∼1-10 kpc, then the wind and the surrounding CGM are both shocked. Warm and cool clouds may be produced in the cooled wind by thermal instabilities, and advected to the halo. Please note that in this scenario the interaction between the hot wind and the ambient medium only provides mass-loading and increases β. The realistic process of the ISM and hot wind interaction needs to be further studied. The central region of the system is zoomed in (see the lower left schematic). The starburst region where the hot wind is launched contributes to the hard X-ray emission, while the interaction between the hot wind and clouds and the large-scale ISM contributes to the soft X-ray emission. If the entrained clouds are cold/cool/warm (∼10 2 -10 4 K), the acceleration of clouds by the hot wind may explain the association between some soft X-rays, Hα or molecular emission. However, it is still unclear whether clouds can be accelerated by the ram pressure of the hot wind to the observed velocities, or they may be completely shredded before fully accelerated due to hydrodynamic instabilities. Magnetic fields may be important to accelerating the entrained clouds. Figure modified from Thompson et al. [186] and Cooper et al. [196].
According to the CC85 model, the momentum deposition rate of a SN-driven wind iṡ Assuming the bolometric luminosity of the host galaxy is L = c 2 SFR with a typical value of ∼ 10 −3 [163,266] for the Salpter initial mass function (IMF), Equation (34) can be rewritten aṡ where −3 = /10 −3 . For the typical value of star-forming galaxies αβ ∼ 1, we haveṀ hot V hot ∼ L/c. In general, the momentum injection by the SN-driven hot wind is comparable to that provided by the radiation field. Therefore, we expect that radiation pressure on dust grains may be important to driving galactic winds. Please note that a radiation-pressure-driven wind is a kind of momentum-driven wind (see the discussion in Murray et al. [251]). We also call these winds as dust-driven winds in order to distinguish them from winds driven by radiation pressure on special lines (i.e., line-driven galactic winds e.g., [267][268][269][270][271][272][273][274]). The relative importance of the SN and radiation feedback have been explored both analytically and numerically. A good example is the numerical simulations by Hopkins et al. [254], who combined SN (including stellar winds and HII) heating and radiation feedback together to model galactic winds in a wide range of galaxies (see Figure 5. They did not include cosmic rays, see Section 4.) Figure 5 shows that the radiation pressure on dust is important in high-redshift and starburst galaxies, while SN heating and momentum injection are dominated in the Milky Way-and Small Magellanic Cloud (SMC)-like galaxies. This is a good demonstration, but more details need to be considered. For example, they used the Tree-SPH code GADGET-III [275], later hydrodynamic simulations of SN-driven winds using sub-grid models for individual SNe have revealed more details of the SN feedback (see Section 2.2). Also, it is assumed that momentum injection by radiation pressure iṡ P rad ≈ (1 + τ IR )L incident , which needs to be studied more carefully by radiation hydrodynamics.
Returning to the theory of stellar winds, the dust-driven winds in the stellar context has been well studied for decades [276][277][278][279][280][281][282][283][284][285][286][287]. It is believed that dust grains can survive and grow around the cool, luminous red supergiants, and AGB stars. Radiation pressure on dust grains can drive stellar winds from these stars. Lamers & Cassinelli [288] reviewed the necessary condition for driving winds by dust. First of all, dust grains need to survive against sublimation. If the radiative equilibrium temperature ≈ (F/a r c) 1/4 with F being the radiation flux onto the dust and a r being the radiation constant is above the "condensation temperature", then the grains will sublimate rather than grow. Different dust grains with different sizes have various sublimation temperature, with a typical value between ∼500-3000 K [289][290][291]. In addition to dust sublimation, dust sputtering also occurs due to high-energy atom/ion incident on a grain surface and could destroy a fraction of dust mass in the environment of hot gas [292][293][294][295][296][297][298][299]. Another necessary condition is that the dust is coupled with the gas. The dust grains obtain momentum from the absorption or scattering of photons, and they collide with the gas [276], then the drag force between the dust grains and the gas causes the coupling of the dust to the gas. The requirement of the momentum coupling sets up the lower limit of the mass loss rate and velocity of a dust-driven wind [300,301]. Similar necessary conditions also need to be considered in the galactic context. Murray et al. [251] estimated that the typical length below which the dust is dynamically coupled with the gas is a few hundred kpc, which is significantly larger than the typical scales on which the galactic winds are launched and accelerated. Although the dust and the gas are dynamically coupled, they may not be in thermal balance. The dust is thermally in equilibrium with the radiation field, but the gas may have different temperature. Goldsmith [302] showed that the dust and gas are only well coupled in regions denser than ∼10 4 -10 5 cm −3 . Please note that the inner few hundred pc of ULIRGs has mean gas density ∼10 3 -10 4 cm −3 . Narayanan et al. [303,304] found that the dust and gas are thermally coupled in ULIRGs. However, the density of the dust-driven winds is much lower than that in the galactic disks, therefore the dust and gas in the winds are dynamically coupled but not thermally coupled (i.e., V dust = V gas but T dust = T gas ). The dust temperature can be estimated by T dust ≈ (F/a r c) 1/4 ≈ 64 K (F/10 13 L kpc 2 ) 1/4 , and the dust is still possible to survive in the hot medium with temperature much higher than the dust sublimation temperature.
Since the dust is coupled with the gas, the dust opacity not only depends on the composition and size of the dust grains, but also depends on the dust-to-gas ratio. Some widely used dust opacities include the frequency-averaged Rosseland and Planck mean opacities [305,306]. Figure 6 shows these two opacities in the literature, most of which was discussed in Semenov et al. [305]. For dust temperature 150 K, the opacity is scaled as κ ≈ κ 0 T 2 . Thompson et al. [252] discussed that the T 2 scaling follows from the fact that the dust absorption cross section scales as λ −2 ∝ T 2 in the Rayleigh limit. The normalization κ 0 depends on the properties of grains and dust-to-gas ratio. For the Galactic environment κ 0 ∼ (2-3) × 10 −4 cm 2 g −1 for the Rosseland mean κ R and κ 0 ∼ 10 −3 cm 2 g −1 for the Planck mean. For dust temperature between 150 K T 1500-2000 K, the dust opacities are almost flat around 5-10 cm 2 g −1 ; and for dust temperature 1500-2000 K, the opacities significantly drop, because the temperature has reached the dust sublimation temperature, but still below the hydrogen ionization temperature ∼10 4 K. Then the gas dominates over the opacities at T 10 4 K.  Figure reproduced using the publicly available code of Semenov et al. [305]. Other opacities [308][309][310][311][312][313][314] were also plotted and compared to Semenov et al.'s opacity in [305]. Please note that the Rossenland mean opacities are scaled as κ ∝ T 2 while the Planck mean opacities are slightly different from κ ∝ T 2 at low temperature. The opacities decrease significantly at the sublimation temperature of dust at T ∼ 1000 K. The rise in the opacity at higher temperature is due first to H scattering, then to bound-and free-free interactions, then to scattering off of free electrons.
Considering gravity, a basic requirement for radiation pressure to drive winds is that the radiation flux of a galaxy near or approach the Eddington limit for dust. A galactic disk which is optical thick to the UV radiation (τ UV > 1) requires the gas surface density Σ g > 1/κ UV , or Σ g > 10 −3 g cm −2 ∼ 5 M pc −2 , where κ UV ∼ 10 3 cm 2 g −1 is a characteristic UV opacity [252]. The typical surface density of rapidly star-forming galaxies and starbursts is above this value and are optical thick to UV. For galaxies which are still optically thin to IR, we can derive the Eddington limit for single-scattering limit, which means that all photons are scattered or absorbed once and then escape the system. Therefore κ ∼ 2/Σ g , and the Eddington limit flux for these galaxies (see details in Andrews & Thompson [307]) is where f g is the mass fraction of gas. For galaxies (mostly starbursts) which are optical thick to IR (Σ g > 5000 M pc −2 f −1 dg, 150 ) with f dg, 150 = f dg × 150 being the dust-to-gas ratio, other two Eddington limit fluxes can be derived, one is for warm starbursts (T dust 150 K) that another is for hot starbursts (T dust 150 K but still below the dust sublimation temperature) that Andrews & Thompson [307] compared the above Eddington limit fluxes to observations and found that normal galaxies are sub-Eddington in general, so the impact of radiation pressure is not significant. Coker et al. [315] found that M82 is sub-Eddington on kpc scales, but the single-scattering Eddington ratio is close to be unity on smaller scales ∼250-500 pc. The starburst cores of Arp 220 may near or reach the Eddington limit for dust [316,317]. Most LIRGs and ULIRGs are most likely to be sub-Eddington to the IR radiation [318][319][320], but they can still be locally super-Eddington due to the density inhomogeneity of the disks [321]. Please note that the estimate about whether galaxies reach the Eddington limit includes uncertainties in the dust-to-gas ratio and the CO/H 2 -HCN conversion factors. Recent numerical simulations also showed that dust-driven winds can be launched even in an initially sub-Eddington system (see Section 3.4). Below we first introduce the analytic models for dust-driven shells in Eddington and sub-Eddington systems in Sections 3.2 and 3.3, then look into the interaction between gas and radiation based on recent radiation hydrodynamic simulations in Section 3.4. The importance of radiation pressure on star cluster scales are discussed in Section 3.5.

Analytical Models for Dusty Shell Acceleration by Radiation Pressure
Ignoring the SN feedback, the dynamics of properties of dust-driven winds can be studied from first principles. The momentum equation of a dust-driven wind is (see more details in [288], Chapter 7) where M is the total mass inside radius r, F is the radiation flux, and κ is the dust opacity. For a spherical wind, the flux can be written as a function of the total luminosity Integrating Equation (41) over dm = 4ρπr 2 dr, and note that the total radiation force is where τ = ∞ R κρdr is the radial optical depth, and R is the starburst region where the dusty wind is launched. As a result, the integrated momentum equation of the wind iṡ whereṀ is the mass-loading rate, which is assumed to be a constant for spherical wind, and v ∞ is the asymptotic velocity of the dusty wind. For Eddington ratio Γ d > 1, the dusty gas is unstable to driving an unbound outflow, but for Γ d < 1, only "fountain flow" can be produced. Ignoring the gravity, we haveṀv ∞ ≈ τL/c. For UV photons, we take the single-scattering limit, i.e., all UV photons only scatter once, thus τ UV ≤ 1. On the other hand, infrared (IR) can be scattered many times and τ IR 1 is possible. Theoretically, the momentum of the dusty wind is boosted from L/c to τ IR L/c L/c for systems which are optical thick to IR. Now we consider a dusty shell accelerated by radiation pressure. Assuming the mass of the shell is M sh , and the shell has a radius of r and thickness of ∆r. The momentum Equation (39) can be rewritten as Here, we take τ = κρ∆r as the optical depth of the shell. If we use an isothermal sphere to model the gas as a function of radius ρ(r) = f g σ 2 /(2πGr 2 ) with f g being the mass fraction of gas and σ being the velocity dispersion, and assume the mass M sh is from the swept-up gas inside radius r, thus, we have M sh as a function of r that M sh = 2 f g σ 2 r/G. Assuming the single-scattering limit τ = 1, we can use Equation (44) to solve the velocity of wind as a function of r: This solution (45) is given by Murray et al. [251]. Here L M is the Eddington luminosity for single-scattering limit that L M = GM sh Mc/r 2 , or in particular, L M = 4 f g cσ 4 /G for an isothermal sphere. An unbound dust-driven wind requires its host galaxy to have a luminosity L > L M , i.e., a super-Eddington luminosity. On the other hand, the L M ∝ σ 4 may be related to the Faber-Jackson relation [322]. Murray et al. [251] argued that a starburst that reaches L M moderates its star formation rate and its luminosity not to increase significantly further. Elliptical galaxies reach L M during their growth at z 1, and that this is the origin of the Faber-Jackson relation.
A more generalized equation for a dusty shell includes both UV and IR opacities. Thompson et al. [323] (see also [11]) modified Equation (44) as For a system with τ UV ≤ 1, the optical depth term 1 + τ IR + e −τ UV goes to τ UV . For τ UV 1 the optical depth term becomes 1 + τ IR . The generalized Eddington limit for a shell is defined as and the Eddington radio Γ = L/L Edd can be reduced to Γ → Γ IR = κ IR L/(4πGMc) for τ IR 1, and Γ → Γ UV = κ UV L/(4πGMc) for τ UV 1. The Eddington limit and Eddington ratio for the single-scattering limit are L single = GMM sh c/R 2 and Γ single = LR 2 /(GMM sh c) respectively. Please note that L single is basically the same as L M in Equation (45), but L M is for an isothermal sphere and L single is for a point source.
For geometrically thin dusty shell with a fixed mass, the opacity of the shell ∝ r −2 and eventually becomes optical thin for UV photons. For Γ IR + Γ single > 1 and Γ UV > 1, the shell can be accelerated first by IR radiation, then by UV radiation. Following Thompson et al. [323], the asymptotic velocity of the shell is where R UV = (κ UV M g /4π) 1/2 is the critical radius where the shell becomes optical thin to the incident UV radiaiton, L 12 = L/10 12 L , κ UV,3 = κ UV /10 3 cm 2 g −1 , and M sh,9 = M sh /10 9 M . There are many simplifications in the picture of shell acceleration. Firstly, the SN feedback is ignored. A more realistic model should include the SN-driven wind combined to the radiation-pressure-driven wind, and the dusty shell may exist in the hot flow. Secondly, the shell is assumed to have a fixed mass. In general, a massive shell sweeps up mass of the ISM and CGM, the mass of the shell increases, and the assumption of a constant M sh breaks down. If the swept-up mass is accumulated in the shell, the shell may be decelerated. Thirdly, the hydrodynamic instabilities are ignored in the analytic models. As discussed in Section 2.4, a hot wind may shred clouds in the ISM, and similar instabilities may occur due to the radiation-shell or shell-ISM/CGM interaction. Also, the radiatively driven instabilities may suppress the coupling between radiation and the dusty gas. Recently, more detailed analytic models of momentum-driven winds have been developed [324]. In addition, multidimensional MHD and radiation hydrodynamic simulations need to be carried out to further explore the dynamics of the momentum-driven shells.

Self-gravitating Luminous Disks
In contrast to the spherical systems, cylindrical disk systems show different properties of dust-driven winds. The Eddington ratio is no longer a constant but varies along the height above the disks. Considering a disk with uniform brightness and surface density I(r ≤ r rad ) = I and Σ(r ≤ r D ) = Σ, where r rad and r D are the outer radius of the luminous and gravitating portion of the disk. Then the gravity along the polar z-axis is (more details see Zhang et al. [325]) (49) and the vertical radiation force is Then the Eddington ratio along the pole is a function of z where Γ 0 is the Eddington ratio at the center of the disk. For r D ≈ r rad , we have Γ ∞ = 2 Γ 0 , and Γ ∞ = 2 if the disk is at Eddington limit. The bright self-gravitating disk is unstable to driving dusty winds. This result is obviously different from that for a spherical system which has a constant Γ and v ∞ ∝ (Γ − 1). For a test particle along the pole, the equation of the particle is wherer = r rad /r D andẑ = z/r D . For Γ 0 = 1 and r D = r rad , the asymptotic terminal velocity of the test particle is where the characteristic velocity of the wind v c = √ 4πGΣr D 500 km s −1 Σ 1/2 0 r 1/2 D,1kpc , Σ 0 = Σ/g cm −2 , and r D,1kpc = r D /1 kpc. Equation (53) gives the maximum asymptotic velocity of a test particle from a disk at the Eddington limit Γ 0 = 1. Moreover, using the Schmidt law Σ SFR ∝ Σ 1.4 gas [326,327] where the gas fraction is f g = 0.5 f g,0.5 . The v ∞ ∝ SFR 0.36 law may be consistent with some observations in ULIRGs and star-forming galaxies [96,104].

Coupling between Radiation Field and Dusty Gas
The optical depth in some dense star cluster and many rapidly star-forming galaxies are thick to the IR radiation, and the dust reprocessed IR radiation may be absorbed and re-emitted multiple times in these systems. Equation (43) shows that the momentum of a dust-driven wind isṀv ∞ ≈ τ IR L/c, where L/c is the initial momentum provided by the radiation field. However, this is only an analytic estimate, in which we assume the wind is spherical and the radiation field is isotropic along each direction. Roth et al. [328] found that the coupling between radiation and gas may be significantly reduced due to the geometry effect: the IR photons may escape along sightline of lower column density of a non-spherical torus, and the momentum transfer from a radiation field to the gas is limited to ∼1-5 L/c even τ IR is big along some special directions.
Importantly, radiatively driven instability e.g., [329,330] may create low-density channels in the gas, and suppress the momentum coupling between radiation and gas because a fraction of the radiation escapes from these channels. It has been argued that the rate of momentum deposition will never exceed a few L/c due to the radiative Rayleigh-Taylor instability (RRTI, [236,331]). To better understand the dynamics of radiation-gas interaction, a series of multidimensional radiation hydrodynamic (RHD) simulations have been carried out. The results depend on the algorithms for solving the RHD equations and may also depend on the spatial resolution of the simulations.
We can simplify a problem as follows: an IR radiation field is injected into a computational box vertically at the bottom of the box and interacting with a shell of dusty gas with a vertically stratified structure, and the gravity is also along the vertical direction. A key question for this problem is: what is the efficiency of momentum transfer from the radiation field to the gas. In other words, the momentum of the gas in the local box is written asṀ what is the efficiency η for a local box simulation? And is it possible to drive an unbound dusty wind by radiation pressure? The results depend on the numerical algorithms. Therefore, we need to briefly look into the methods to solve the RHD equations. The time-dependent radiative transfer equation is where ν is the frequency, I ν is the radiation intensity, η ν is the emissivity, and κ ν is the opacity.
Integrating Equation (56) over a solid angle and the frequency, we obtain the zeroth and first momentum equations where are the frequency-integrated gray energy, flux, and tensor of the radiation field. Also, the frequency-integrated energy and momentum source term S r (E) and S r (P) can be added to the hydrodynamic energy-and momentum-conservation equations respectively (see [332,333] for more details. The hydrodynamic equations coupled with radiation, combined with the radiative transfer Equation (56) or the zeroth and first momentum Equations (57) and (58) give the set of RHD equations.
Since the total number of unknown fluid and radiation variables are more than the number of equations, we need other inputs, i.e., the so-called closure relations relating P to E r to solve all the variables [334]. A widely used closure relation is the flux-limited-diffusion (FLD) approximation. Neglecting the first momentum Equation (58), and assume the radiation field is locally isotropic, so P = E r I/3, and where λ(E) is the flux-limiter to prevent the flux becoming unphysical [335]. As a consequence, the FLD approximation (62) with the zeroth momentum equation closes the RHD equations. A better treatment than the FLD approximation is the M1 closure method. The Eddington tensor is defined as f ≡ P r /E r . The key technique is to calculate the Eddington tensor and close the set of RHD with both the zeroth and first radiation momentum equations. The M1 closure method calculates the Eddington tensor as where n = F r /|F r |, also χ depends only on the reduced flux and can be calculated following [336]. A more elaborate but more expensive closure method is called the variable Eddington tensor (VET) method, in which the Eddington tensor is direct calculated by the definition in each timestep: and the intensity I is solved based on the time-independent radiative transfer equation using the short characteristic method [337][338][339][340][341].
Krumholz & Thompson [344] used the FLD algorithm to investigate the efficiency of momentum coupling between IR radiation and vertically stratified gas in a local 2D box and found that the radiation pressure on dust drives RRTI and supersonic turbulence. Although the radiation trapping factor η can be big at the base of the system, a system which is initially sub-Eddington for dust can never launch an unbound wind. In the following work, Krumholz & Thompson [345] turned off the gravity and prolong the 2D computational box to investigate the maximum efficiency of radiation-gas coupling using the FLD algorithm. Without gravity, a wind can be launched and accelerated, but the RRTI-like instability (without gravity, only drives by the radiation pressure) forces the gas to spread out and forms many low-density channels, which limit the momentum transfer from the radiation to a few L/c, regardless of the optical depth of the system. Thus, the efficiency η is low, and lower for higher τ IR . They concluded that an initially sub-Eddington system can never launch and accelerate winds-since most ULIRGs are below the IR Eddington limit for dust, radiation pressure is not likely to be crucial to galactic winds.
The problems studied by Krumholz & Thompson [344] (hereafter the KT problems) have been revisited by many others using different algorithms based on various RHD codes. Rosdahl & Teyssier [347] adopted the M1 closure method to study the KT problems. They found that the dusty gas receives a larger acceleration than in the FLD simulations, but the acceleration is not sufficient to overcome the gravity, and the gas eventually settles down at the bottom of the system, similar to the results in Krumholz & Thompson [344]. Kannan et al. [349] also used the M1 closure method to re-investigate the KT problems and found similar results.  middle and right panels). The initial IR optical depth for dust is τ IR = 3 (left and middle panels), and τ IR = 1 (right panel), with the same radiation field injected at the bottom boundaries. The FLD run shows that radiatively driven instability creates low-density channels and the dusty gas spreads out of the entire box. The VET runs show that the dusty gas is still accelerated as a whole with more IR flux trapped in the gas, and more efficient momentum coupling between radiation and gas than that in the FLD run.
Using the VET algorithm, Davis et al. [346] also revisited the KT problems. In contrast to the results using the FLD and M1 closure method, Davis et al. showed a stronger momentum coupling between radiation and dusty gas. Importantly, although the RRTI develops and limits the radiation-gas interaction, an initial sub-Eddington system can turn to a super-Eddington system due to the radiation heating on dust, and the gas can still be accelerated upward by radiation and produce an unbound outflow. In the absence of gravity, Zhang et al. [319] using the VET method showed that the gas spreads out along the vertical direction, but the gas is still accelerated as a whole. In contrast to the result in Krumholz & Thompson [345], Zhang et al. found that the momentum transfer from the radiation to the gas is not merely a few L/c, but the efficiency η ∼ 0.5-0.9 for a large range of τ IR . Therefore, the moderate resultṀv ∞ ∼ [1 + (0.5 − 0.9)τ IR ]L/c is between the analytic estimate for spherical wind η = 1 and the FLD result η 1. Figure 7 shows some detailed comparison between a FLD run and two VET runs. In summary, the RHD simulations using the VET algorithm showed more efficient momentum boost than those using the FLD and M1 closure method. Also, similar to the VET simulations, the RHD simulations based on the Monte Carlo radiation transfer scheme shows similar results [348]. Different from Krumholz & Thompson [345], Zhang et al. proposed that radiation pressure is still important to driving dusty winds in most rapidly star-forming galaxies which are initially sub-Eddington but optically thick to the IR radiation.
In addition to the algorithm-dependent differences, the behavior of the radiation-pressure-driven dusty shells also depends on the spatial resolution of the simulations. The KT problems focus on the resolution of the pressure scale height with a resolved length to h * = c 2 s /g with c s being the sound speed in the gas, and the resolution ∼h * is far below the pc scales. Recently, Krumholz [350] studied the direct stellar radiation interacting with the dusty accretion flows and found that whether radiation feedback is able to reverse an accreting/falling flow depends on the spatial resolution of the simulations. They argued that the IR radiation pressure on dust may also be overestimated by low-resolution simulations. On the other hand, it has been tested that the behavior of a dusty shell on the scales of pressure height h * is very different from that on galaxy scales with lower resolution. The study on the effects of spatial resolution is still underway.

Radiation Feedback in Star Clusters
Even galaxies are below the average Eddington limit for dust, their star clusters which are much denser than the surrounding medium may still reach or exceed the Eddington limit on small scales, e.g., [351]. The models and calculations in Sections 3.2-3.5 can be applied to star clusters. Murray et al. [253] proposed that the radiation pressure from clusters with mass 10 6 M is able to expel the surrounding dusty gas out of the clusters and the galactic disk with velocities ∼ a few× 100 km s −1 . Thus, the cold/cool outflowing gas can be accelerated by radiation from the disk, and travel to a distance of ∼50-100 kpc. These outflows may be observed as the Mg II or Na D absorbers. Zhang et al. [320] carried out RHD simulations to model the dynamics of dusty clouds accelerated by a radiation field, and found that the dusty clouds can be accelerated to hundreds of km s −1 with a significant longer lifetime compared to those entrained in a hot flow where the momentum injection is comparable to that in the radiation field.
In addition to driving galactic winds, radiation pressure from star clusters also disrupts the parent GMCs, drives turbulence and affects star formation before the explosion of the first SN. In particular, multidimensional RHD simulations have been carried out to study this problem. Skinner & Ostriker [258] performed simulations using the M1 closure method: the IR flux is injected from super star clusters with τ IR > 1 and interacts with a turbulence GMC. They found that radiation pressure reduces the star formation efficiency only if the dust opacity κ IR 15 cm −2 g. Tsang & Milosavljević [262] did similar simulations using the Monte Carlo radiation code and basically reaffirm the main conclusions in Skinner & Ostriker [258]. Raskutti et al. [259,261] also used the M1 closure method and found that the direct FUV radiation pressure from star clusters with τ IR < 1 is also important to suppressing the star formation efficiency in GMCs.

Some Early Work and Analytic Work
Stellar feedback also includes cosmic rays (CRs). The shocks of SN remnants are efficient accelerators of CRs by diffusive shock acceleration (first-order Fermi acceleration) [352][353][354][355][356][357][358][359]. About ∼10% of the SN energy (∼10 50 ergs) can be converted into non-thermal CR energy [360][361][362], most of which is deposited in protons with their energies following a power-law distribution peaked at ∼ GeV [363][364][365][366]. The local energy density of CRs is ∼1.8 eV cm −3 [367], and the total energy of CRs in our Galaxy is ∼10 56 ergs. It appears that the GeV CRs are confined to the Galaxy for about 2 × 10 7 yr [368,369], therefore the CR luminosity is ∼10 41 ergs s −1 in our Galaxy. It has long been proposed that CRs coupled to the gas may also drive large-scale galactic winds in the Galaxy.
Ipavich [370] first studied the CR-driven galactic winds. He found that if CRs travel through a magnetized plasma medium faster than the Alfvén speed, they can be coupled to the plasma by the emission of MHD waves. The mass-loading rate of CR-driven winds can reach ∼1-10 M yr −1 . Breitschwerdt et al. [371] studied the interaction of CRs and hot gas along large-scale magnetic fields and found that the combination of CR and thermal effects can drive a supersonic galactic wind with a mass loss rate of ∼1 M yr −1 in the Galaxy. Breitschwerdt et al.'s work was further extended by their group [372][373][374][375]. For example, Zirakashvili et al. [373] suggested that the coupling between CRs and plasma is provided by the resonant excitation of small-scale magnetic field fluctuations generated by CR streaming instability [376,377], and the nonlinear Alfvén wave damping is sufficiently strong to heat the plasma. They investigated the winds from a rotating disk galaxy and found that the mass loss rate is larger compared to that from a non-rotating and undamped model. The winds also take away a large fraction of angular momentum. Ptuskin et al. [374] further explored CR streaming instability and the transport of relativistic nucleons in CR-driven winds. Everett et al. [378,379] considered galactic winds driven by the combination of both CRs and thermal gas pressure under Milky Way conditions and found that both CRs and thermal gas pressure are essential to explaining the observed Galactic diffuse soft X-ray and radio emission.
For star-forming and starburst galaxies, the total energy injected by CR protons is [380,381] L CR ≈ 3.2 × 10 41 ergs s −1 e 51 µ 100 SFR 1 ≈ 6 × 10 −4 e 51 µ 100 where e and µ have been given in Section 2.2, and L is the galactic luminosity. Although L CR L, the momentum injection by CRs may still be significant. CRs scatter many times due to the interstellar magnetic irregularities. The CR escaping timescale in the Galaxy is about 2 × 10 7 yr, thus the mean free path of the GeV CRs is λ CR ∼ 1 pc, corresponding to a diffusivity of 10 28 -10 29 cm 2 s −1 . Some other work suggested that λ CR ∼ 0.1 pc [382]. For a typical scale height of winds H ∼ 1 kpc, the CR scattering optical depth is τ CR ∼ H/λ CR ∼ 10 3 − 10 4 , and the momentum deposited by CRs is which is comparable to the momentum injected by the direction radiation. Socrates et al. [381] argued that the CR feedback may be important in star-forming galaxies and starbursts to driving galactic winds. However, this estimate includes several uncertainties. The mechanism responsible for the mean free path λ CR and CR diffusivity in our Galaxy are not fully understood, and λ CR in other galaxies are not known. CRs are destroyed by interacting with ambient protons and producing pions with a timescale of t pp ∼ 7 × 10 7 yr/(n H / cm −3 ) [380,383]. For starbursts with n H ∼ 10 3 − 10 4 cm −3 , t pp ∼ 7 × 10 3 − 10 4 yr may be significantly shorter than the CR diffusion time, therefore the hadronic losses cannot be ignored. More recently, it has been found that the momentum injection from an individual SNR to the ISM can significant increase once the effect of CRs produced by the SNR is also included [384]. Numerical simulations to explore stellar feedback with the combination of SNRs and CRs are underway. In addition, how to model the impact of CRs in star-forming galaxies and starbursts is still an open question.

Recent Numerical Simulations
Recently, 3D hydrodynamic and MHD simulations have been performed to study the CR-driven galactic winds, with various approximations on CR transport [385][386][387][388][389][390][391][392][393][394][395][396][397][398][399][400][401][402]. Since the properties of the CR-driven winds sensitively depend on the details of CR transport, we first look to the set of two-fluid MHD equations which describe the winds where ρ, u, B and g are the gas density, velocity, the local magnetic field and gravity, respectively. The total pressure, MHD pressure and the total MHD density are with P th , P cr being the thermal and CR pressure, and P A is contributed by the Alfvén wave. The CR energy source term Q cr is the heating rate from SNe (∼10% of the SN energy), Q had is the hadronic losses [393,398], and Q Coul is the Coulomb losses due to the CRs interaction with the ambient gas [403].
Other variables are discussed as follows: (1) Other Sources. The mass, momentum and energy injection rates by other sources depend on the models of galaxies. Hereṁ other is the combination of the mass injection rate by sources including jets, stellar and SN-driven winds subtracting the mass used for star formation e.g., [400]. The momentum injection rateṖ other can be provided by SNe or radiation in the absence of AGNs. The energy rateu other is also given by SNe, and radiative cooling needs to be subtracted from this term. (2) Alfvén Wave Damping. As discussed by early theoretical work e.g., [371,373,374,378,379], the energy source term S A is caused by Alfvén wave damping, which is still not completely understood [366,382,401]. The damping mechanisms include the ion-neutral damping, which is caused by the friction between ions and neutrals [404][405][406][407], linear and nonlinear Landau damping [408][409][410][411][412][413], and turbulent damping [365,366,[413][414][415][416]. Everett & Zweibel [417] considered both ion-neutral and nonlinear Landau damping, and found that these damping mechanisms are only important if the magnetic fields are above 10 µG, or high CR pressure (∼10 −11 ergs cm −3 ). The CR diffusion also depends on the wave damping.
CR Diffusion. The scattering term ∂e cr /∂t| scatt in Equation (72) includes CR diffusion and streaming. The diffusion term is usually written as ∇ · (κ · ∇e cr ), where κ is the diffusion coefficient. For isotropic CR diffusion, the diffusion term can be rewritten as [418] ∇ · (κ · ∇e cr ) = ∇ · (κ bb · e cr ) + ∇ · [κ ⊥ (I − bb) · e cr ]. (76) where κ is the diffusion tensor, the κ and κ ⊥ are the diffusion coefficients parallel and perpendicular to the magnetic field, and b = B/|B| is the unity vector along the magnetic field.
If magnetic field lines are sufficiently tangled on small scales, CR diffusion can be approximated to be isotropic [418]. As mentioned in Section 4.1, the value of diffusion coefficient in our Galaxy can be estimated by κ ∼ H 2 /t CR with H ∼ 1 kpc and the CR traveling time t ∼ 10 7 yr so κ ∼ 3 × 10 28 cm 2 g −1 . However, CR diffusion is considered to be anisotropic in most cases, and the coefficient κ ⊥ is found to be much lower than the parallel one κ κ ⊥ . The diffusivity also depends on the rigidity and may vary in the multiphase medium. The rigidity is R i = r L Bc with r L being the gyroradius, and κ ∝ R 0.6 i [419,420]. Farber et al. [400] suggested that κ depends on the gas temperature and used κ = 10 29 cm 2 s −1 for T < 10 4 K, and κ = 3 × 10 27 cm 2 s −1 for T 10 4 K.
On the other hand, it is unclear if the Galactic diffusivity can be applied to other galaxies. One example is NGC 253. Heesen et al. [421][422][423] investigated the magnetic fields in NGC 253 and estimated that the average diffusion coefficient is κ ∼ 2 × 10 29 cm 2 s −1 , and κ ⊥ ≈ 1.5 × 10 28 cm 2 s −1 E(GeV) 0.5±0.7 , which seems to be one order of magnitude higher than the values in our Galaxy. The comparison between the observed values for NGC 253 and the theoretical model in [424] was discussed in [425]. The estimate of the CR diffusivity in other galaxies has also been explored [426].
A widely used CR streaming model for galactic wind simulations calculates the streaming velocity as v st = − − v A (B · ∇P cr )/|B · ∇P cr | with v A being the Alfvén velocity, and the energy term ∂e cr /∂t| scatt contributed by streaming is −|v A · ∇P cr |. The CR energy flux due to streaming is given by F cr = v st (E cr + P cr ), and the total CR energy flux including the diffusive flux along the direction of local magnetic field is F cr = (E cr + P cr )(v + v st ) − bκ(b · ∇P cr ). In the self-confinement scenario, in which CRs scatter on waves excited by the stream instability [365,376,377,382], an effective draft speed f v st is used to replace v st , and the factor f can be calculated by balancing the wave growth rate with the wave damping rate [394,427,428].
Another way to calculate F cr is to directly solve it by adding a new equation and some closure relations to the set of Equations (68) to (72). Before we introduce the more sophisticated algorithms to solve the equations for CR-driven winds, we briefly review a series of recent 3D global and local stratified numerical simulations in the literature. All the simulations include some physical processes introduced above but ignored others. Most simulations showed positive results that the CR feedback drives galactic winds in various galaxies. Uhlig et al. [385] studied the impact of CR streaming at the local sound speed, and found that powerful winds are observed in dwarf galaxies (∼10 9 M ) and low-mass galaxies ( 10 11 M ), while fountain flows are driven in larger halos ( 10 11 M ). This halo-dependent conclusion is confirmed by Jacobi et al. [396]. Hanasz et al. [387] considered anisotropic diffusion of CRs in star-forming (40 M yr −1 ) disk galaxies with gas surface density Σ g ∼ 100 M pc −2 , and found that CRs alone can trigger the formation of strong winds. Booth et al. [386] assumed isotropic diffusion and compare their simulations to Milky Way and Small Magellanic Cloud (SMC)-like galaxies, and found that the mass-loading rate is up to β CR ∼ 10 in dwarf systems, in which warm gas ∼10 4 K is dominated the winds. Salem & Bryan [388] also assumed isotropic diffusion and found stable winds with β CR ∼ 0.3 from ∼10 12 M halos. Interestingly, Fujita et al. [397] did similar simulations as in Salem & Bryan [388] but found no significant difference in mass-loading rates of SN-driven winds with or without CR pressure. Pakmor et al. [390,391] argued that anisotropic diffusion is more reliable than the isotropic approximation for galaxies with 10 11 M halos. Girichidis et al. [389] performed the first MHD simulations in stratified boxes to model CR-driven winds from the Milky-Way-like disks Σ g ∼ 10 M pc −2 with SNe, assuming anisotropic diffusion for CRs. They also included hadronic cooling in [398] and found that 5-25% of the injected CR energy cools via hadronic losses. Ruszkowski et al. [394] included both the CR (self-confinement) streaming and anisotropic diffusion in their MHD model and found that the presence of moderately super-Alfvén CR streaming enhances the efficiency of galactic winds from galaxies with 10 12 M halos. Farber et al. [400] extended Ruszkowski et al. [394]'s work by considering temperature-dependent anisotropic diffusion and observed the decoupling of CRs in the cold/cool neutral ISM (<10 4 K). CRs propagate faster in the cold/cool gas than in the ionized medium. Holguin et al. [402] further developed a more elaborate model on CR streaming u st = f u A with f depending on the turbulent Mach number, and investigated the effects of turbulence damping. Wiener et al. [395] investigated the relative importance of CR diffusion vs streaming and found significant differences between these two mechanisms in dwarf galaxies (∼10 10 M ). On the other hand, the heating by CR-driven winds in CGM and clusters of galaxies has also been explored e.g., [403,[429][430][431] We expect that future numerical simulations of the CR feedback may use more sophisticated algorithms for CR transport. Returning to the CR Vlasov equation, the advection-diffusion equation for CR distribution function is [399,432] Equation (77) (77) is the CR energy Equation (72), and the first momentum equation using the reduced speed of light approximation is Please note that we have reduced their P cr tensor to P cr , κ is the "equivalent" diffusion coefficient, and the technique to use the reduced speed of light (c < c) can be also found in RT algorithms [320,342,433]. As a result, the CR streaming energy flux no longer holds F cr = (E cr + P cr )(v + v st ), but Equation (78) is added to the original set of MHD equations to solve the CR transport. Jiang & Oh [399] found that their two-moment method is more stable and robust to handle CR streaming. Recently, Thomas & Pfrommer [401] used the Eddington approximation of RT to develop another two-moment algorithm of CT transport in the self-confinement picture, including both the CR streaming, diffusion, and other damping mechanisms. The CR energy density, flux and Alfvén wave density can be solved together by the set of MHD equations. Compared to Jiang & Oh [399]'s method which is correct to the order of v A /c, the CR transport algorithm by Thomas & Pfrommer [401] reaches the order of (v A /c) 2 . These two-moment algorithms are potentially better to model the CR feedback and CR-driven winds.
One may ask whether CRs can accelerate cool/cloud clouds [417,428,434]. It is found that the CR pressure gradient can push the cool/cloud clouds outward, but CRs can be decoupled with cool/cold gas. It seems that the velocities gained by the pressure gradient is not sufficient to match the observed velocities [428]. Please note that most of the CR-driven wind simulations are on galactic scales, it is also important to perform small-scale simulations to better understand the interaction between CRs and multiphase gas especially cool/cold clouds.

Conclusions and Discussion
I have reviewed the theory of galactic winds driven by stellar feedback processes, including supernova (SN) explosions, radiation pressure from starlight on dust grains, and cosmic rays.
Core-collapse SN explosions have long been considered as the primary energy and momentum sources of galactic winds. A SN-driven large-scale galactic winds can be approximated described by the 1D spherical solution of the CC85 model [161], which is controlled by two key parameters: the thermalization efficiency that measures the thermal energy converted from the net SN explosions to the wind, and the mass-loading rate. These two parameters can be constrained using the X-ray data from star-forming and starburst galaxies, or by hydrodynamic simulations. Please note that the CC85 model is for adiabatic flows, radiative cooling may significantly change the properties of winds above a distance of ∼1-10 kpc.
A prevailing scenario of the multiphase galactic winds is that the multiphase material is advected into a SN-driven hot wind and accelerated by the ram pressure of the hot wind to a velocity of a few hundred km s −1 . This scenario is supported by some global simulations on galactic scales and some observations of the association between multiphase gas, but recent numerical simulations on smaller scales found that clouds entrained in a hot flow are completely shredded via hydrodynamic instabilities long before being fully accelerated by the hot wind. Magnetic fields may be important to prolonging the lifetime of clouds, but more work needs to be done to explore the cloud dynamics in hot flows as well as large-scale interaction between hot flows and the multiphase ISM in a wide range of galaxies with more realistic structure of magnetic fields. Another possible origin of warm/cool outflowing gas arises from thermal instability of the hot winds.
On the other hand, radiation pressure from starlight on dust grains can be important in rapidly star-forming and starburst galaxies. Massive stars emit UV photons, which are absorbed by dust and re-radiated in IR from dust grains. Many galaxies are above the Eddington limit for UV photons, but only galaxies with flux 10 13 L kpc −2 may reach the Eddington limit for IR photons. The estimate of the Eddington limit is uncertain in the dust-to-gas ratio. However, even a galaxy is globally sub-Eddington for dust, its star clusters which are much denser than the surrounding medium may still reach the Eddington limit, and drive dusty outflows from them. The dusty gas from star clusters may be eventually pushed out of the galaxies by the momentum injection of SNe and radiation pressure. The analytic models have been proposed to describe the geometrically thin dusty shells pushed by radiation pressure, while the behavior of dusty gas from a self-gravitating disk is different from that from a spherical system. The momentum coupling between the dust and UV photons can be modeled using the single-scattering limit, in which each photon is only scattered or absorbed once. The coupling between the dust and IR photons can be given by P rad ≈ (1 + ητ IR )L/c, with η ∼0.5-0.9, obtained from recent radiation hydrodynamic simulations using the state-of-the-art algorithms. According to the simulations, an initially sub-Eddington system may still launch an unbound wind by radiation pressure on dust. In addition to driving galactic winds, radiation pressure can also disrupt GMCs, which is beyond the scope of this review article.
Whether CRs are important to driving winds in star-forming and starburst galaxies is less unclear, as most of the CR knowledge is gained from our Galaxy. For example, the CR diffusivity may be different between our Galaxy and starbursts such as NGC 253. Taking the CR properties in Milky Way as the fiducial values, CRs are no less important than SNe or radiation pressure to providing momentum injection and coupling to the gas in various galaxies. Cool/cold gas may be decoupled with CRs but still be accelerated by the pressure gradient generated by CR transport. Recently, a series of hydrodynamic and MHD simulations have been carried out to explore the impact of CR streaming and diffusion, primarily focusing on Milky Way-like galaxies, then extending to star-forming and high-redshift galaxies. Two-moment algorithms used for radiative transfer may also be applied to solve the CR transport equations. We expect the two-moment algorithms can be used to model the CR-driven winds in the near future. Also, the impact of CRs on rapidly star-forming and starburst galaxies need to be explored more detailed. A more realistic model for galaxy evolution and galactic winds needs to combine all feedback mechanisms including stellar and black hole/AGN feedback together, as many galaxies include both rapidly star-forming regions as well as central AGNs. For stellar feedback alone, some work has already combined SN and radiation pressure feedback, or SN and CR feedback together, and attempted to build a unified model for stellar-driven galactic winds. However, the details of galactic winds need to be investigated systematically and compared to the multiwavelength observations more carefully.
For example, can the momentum flux from a momentum-driven wind be boosted by an optically thick system with multiple scattering? Although we found that numerical algorithm is important to model the dynamics and thermal properties of momentum-driven winds, recently new evidence showed that spatial resolution may be crucial to modeling SN and radiation pressure, and sub-grid models need to be adopted for numerical simulations on galactic scales. Also, geometry effect may be important: feedback simulation in a vertically stratified box shows different result from that in a global simulation with the same setup. It has been found that the momentum injection by stellar feedback may drive turbulence as well as launching and accelerating winds, so how much momentum is transferred to turbulence and how much to galactic winds? Among SNe, radiation feedback and gravity, which is the most important mechanism to drive turbulence and regulate star formation on various scales? These questions need to be further studied in the near future.
Funding: This research was supported by NSF grant AST-1616171.