Hydrodynamic Simulations and Time-dependent Photoionization Modeling of Starburst-driven Superwinds

Thermal energies deposited by OB stellar clusters in starburst galaxies lead to the formation of galactic superwinds. Multi-wavelength observations of starburst-driven superwinds pointed at complex thermal and ionization structures which cannot adequately be explained by simple adiabatic assumptions. In this study, we perform hydrodynamic simulations of a fluid model coupled to radiative cooling functions, and generate time-dependent non-equilibrium photoionization models to predict physical conditions and ionization structures of superwinds using the MAIHEM atomic and cooling package built on the program FLASH. Time-dependent ionization states and physical conditions produced by our simulations are used to calculate the emission lines of superwinds for various parameters, which allow us to explore implications of non-equilibrium ionization for starburst regions with potential radiative cooling.


Introduction
Thermal and mechanical feedback from OB stars in stellar clusters displaces the surrounding medium in starburst regions on a large scale and forms a galactic-scale outflow named superwind (Heckman et al. 1990), accompanied by a narrow shell and sometime by a hot bubble called superbubble (Weaver et al. 1977). The physical properties of the expanding wind region prior to the bubble and shell have been obtained by Chevalier & Clegg (1985) using adiabatic fluid equations that yield the outflow density n ∝ r −2 and temperature T ∝ r −4/3 . However, the fluid equations coupled to radiative cooling functions studied by Silich et al. (2004) depict a deviation from the adiabatic temperature, which may explain strong cooling and suppressed superwinds seen in observations of some star-forming galaxies Turner et al. 2017;Jaskot et al. 2017). In particular, semi-analytic studies and hydrodynamic simulations demonstrated that radiative cooling is heavily dependent on the metallicity, mass-loss rate, and wind velocity (Silich et al. 2004;Tenorio-Tagle et al. 2005;Gray et al. 2019a;Danehkar et al. 2021). In the case of starburst galaxies where metallicity is low, high mass-loss rates and low outflow velocities contribute to substantial radiative cooling (Danehkar et al. 2021).
Photoionization calculations were performed to identify the superwind models with strong radiative cooling (Gray et al. 2019a;Danehkar et al. 2021). However, emission lines in photoionization (PIE) and collisional ionization equilibrium (CIE) calculated by Danehkar et al. (2021) did not make a clear distinction between those with and without substantial radiative cooling. Previously, photoionization models built with timedependent non-equilibrium ionization (NEI) states by Gray et al. (2019a) and Gray et al. (2019b) also indicate that ions such as O vi and C iv could behave differently where plasma is in the NEI case. Non-equilibrium conditions occur when the radiative cooling timescale τ cool is shorter than the CIE timescale τ CIE (see e.g., Gnat & Sternberg 2007)). In the expanding wind region where plasma is in transition from CIE to PIE at temperatures below 10 6 K, NEI conditions may emerge (Vasiliev 2011) and NEI states could have substantial deviations from CIE states (Gnat & Sternberg 2007;Vasiliev 2011;Oppenheimer & Schaye 2013), which can be identified using the C iv, N v, and O vi lines.
Recently, Danehkar et al. (2021) reported emission line fluxes calculated based on CIE+PIE assumptions using the physical conditions obtained from hydrodynamic simulations with the maihem module in the flash program. Similarly, we also computed emission line fluxes following Gray et al. (2019a) for NEI conditions using the timedependent NEI states and physical properties predicted by our hydrodynamic simulations (Danehkar et al. 2022), showing that the C iv and O vi lines have different behaviors in non-equilibrium photoionized expanding wind regions.

Hydrodynamic Simulations
To study starburst-driven superwinds, we consider starburst feedback from a spherically symmetric stellar cluster parameterized by the cluster radius R sc , the mass-loss rateṀ , and the stellar wind velocity V ∞ . The radiation field is characterized by the total stellar luminosity and spectral energy distribution (SED). The surrounding medium has a density n amb , while its temperature T amb is dependent on the radiation field and determined by a cloudy model.
Our hydrodynamic simulations are performed using the maihem atomic and cooling package (Gray et al. 2015;Gray & Scannapieco 2016;Gray et al. 2019b) in the frame work of the flash program (Fryxell et al. 2000), which obtains the solutions for the following one-dimensional spherically symmetric fluid equations coupled to the radiative cooling and photo-heating functions: where r is the radius, ρ the density, u the velocity, P the thermal pressure, E the total energy per unit mass, γ = 5/3 the specific heat ratio, q m =Ṁ /( 4 3 πR 3 sc ) and q e = ( 1 2Ṁ V 2 ∞ )/( 4 3 πR 3 sc ) the mass and energy deposition rate per unit volume, respectively, n i the number densities of ions, n e the electron number density, Λ i the radiative cooling rates for a specified temperature from Gnat & Ferland (2012) dν the photo-heating rates obtained from the given radiation field J ν and the photoionization cross-section σ i (ν) (Verner & Yakovlev 1995;Verner et al. 1996), ν the frequency, ν 0,i the ionization frequency, and h the Planck constant.
To set the boundary conditions, we employ the semi-analytic radiative assumptions adopted by Silich et al. (2004), which are based on the adiabatic solutions obtained by Chevalier & Clegg (1985). Accordingly, the density, temperature, and velocity at r = R sc  The radiation field J ν included in our hydrodynamic simulations is generated by the stellar population synthesis program Starburst99 Leitherer et al. 2014) for the rotational stellar models Georgy et al. 2012) with an initial mass function with slope α = 2.35 within the range 0.5-150 M ⊙ and the total stellar mass M ⋆ = 2 × 10 6 M ⊙ , which are associated with the mass-loss ratė M = 10 −2 M ⊙ yr −1 at 1 Myr. The Starburst99 radiation field is applied to the photoheating rates in maihem to perform non-equilibrium calculations, as well as cloudy photoionization models. Figure 1 (left panels) presents the temperature T and density n radial profiles (solid red lines) generated by our maihem simulation for a model with substantial radiative cooling, along with the expected adiabatic solutions without radiative cooling (dashed lines). The four distinctive regions of a typical superwind defined by Weaver et al. (1977) are also separated by dotted, dashed, and dash-dotted vertical gray color lines, namely expanding wind region (before dotted), bubble (between dotted and dashed), shell (between dashed and dash-dotted), and ambient medium (after dash-dotted vertical lines). The Strömgren sphere (solid vertical gray color line) are determined by a cloudy photoionization run on the density profile (pure PIE) following Danehkar et al. (2021). Figure 2 shows the mean radiative temperature over the mean adiabatic temperature, f T ≡ T wind /T adi , of the expanding wind region predicted by our maihem hydrodynamic simulations for different wind parameters (V ∞ andṀ ), ambient densities (n amb ), metallicity (Z/Z ⊙ ), a stellar cluster with R sc = 1 pc and M ⋆ = 2 × 10 6 M ⊙ , and current age t = 1 Myr. The catastrophic cooling (CC) and catastrophic cooling bubble (CB) wind modes, which are with and without bubbles, have f T < 0.75, while the adiabatic bubble (AB) and pressure-confined (AP) mode have 0.75 < f T < 1.25. Moreover, the adiabatic Figure 2. The mean radiative temperature T rad with respect to the mean adiabatic temperature T adi of the expanding wind region for various wind parameters V∞ = 500, and 1000 km s −1 anḋ M = 10 −2 ×(Z/Z⊙) 0.72 M⊙ yr −1 , metallicity Z/Z⊙ = 0.125, 0.25, 0.5, and 1, and ambient media with log n amb = 0, 1, 2, and 3 cm −3 , surrounding a stellar cluster characterized by Rsc = 1 pc, M⋆ = 2 × 10 6 M⊙, and age t = 1 Myr. The wind models are identified as adiabatic bubble (AB), catastrophic cooling (CC), catastrophic cooling bubble (CB), and pressure-confined (AP), with optically-thick status (bold font), based on criteria defined by Danehkar et al. (2021). pressure-confined (AP) mode is assigned to those models where the bubble expansion is confined by the ambient thermal pressure (for more detail see Danehkar et al. 2021). Optically-thick models having neutral ambient medium are also displayed with the bold font. It can be seen that increasing mass-loss rates and decreasing wind velocities result in enhanced radiative cooling in the expanding wind region.

Time-dependent Photoionization Models
The physical conditions and NEI states of the expanding wind region made by our hydrodynamic simulations with maihem, along with the radiation field generated by Star-burst99, are used to construct non-equilibrium photoionization models with the cloudy program (Ferland et al. 2013(Ferland et al. , 2017. Excluding the NEI states, Danehkar et al. (2021) incorporated the physical properties into a grid of cloudy models, which describe the CIE+PIE cases. Inclusion of the NEI states predicted by maihem allows us to emulate non-equilibrium photoionization in the expanding wind region (Danehkar et al. 2022) for which pure PIE is still present in the ambient medium.
In the non-equilibrium conditions, gas kinematic and ionization structures are closely interconnected. The time-dependent NEI states are linked to collisional and dielectronic recombination, collisional ionization, and photoionization rates. The number density of ions n i for each chemical element in the non-equilibrium cases can be described as where α i includes the radiative recombination rate (Badnell 2006) and dielectronic recombination rates (references given in Table 1 of Gray et al. 2015) for the ionic species i, S i is the collisional ionization rates (Voronov 1997 The non-equilibrium cases appear when the CIE timescale τ CIE ≈ 1/(n e α i + n e S i ) (Mewe 1999) is longer than the cooling timescale τ cool = 3(n i + n e )k B T /(2n 2 i Λ i ) (Dopita & Sutherland 2003). For less dense environment ( 1 cm −3 ) that is typical of the expanding wind region, the ions C iv, N v, and O vi satisfy the condition τ CIE τ cool at temperatures below 10 6 K, so they could be in the NEI situations in the presence of strong radiative cooling. For τ CIE ≪ τ cool , plasma is in the CIE conditions. Figure 1 (right panels) presents the emissivities of Hβ, low-excitation [S ii] and [N ii], and high-excitation He ii, [O iii], and C iii], as well as highly-ionized C iv and O vi calculated by cloudy for the physical properties (without the NEI states) produced by our maihem simulation associated with plasma in photoionization and collisional ionization equilibrium (top panel; CIE+PIE), as well as the emissivities computed using cloudy following the method of Gray et al. (2019a) for the physical conditions and NEI states generated by our maihem simulation corresponding to the non-equilibrium photoionization situations (bottom panel; NEI with pure PIE in the ambient medium). Large grids for various model parameters are provided as interactive figures by Danehkar et al. (2021) for combined CIE+PIE conditions and Danehkar et al. (2022) for combined NEI+PIE situations, and hosted on this website https://galacticwinds.github.io/superwinds/. As seen in Figure 1, the O vi and C iv emissivity profiles in NEI are not the same as those with CIE (see green and orange dashed lines), particularly in the expanding wind region affected by radiative cooling. This behavior can be explained by the time-dependent ionization states of the ions O vi and C iv at temperatures below 10 6 K when rapid cooling faster than ionization processes occurs (τ cool < τ CIE ).

Implications for Starburst Galaxies
Our NEI calculations indicate that radiative cooling could enhance the C iv 1550Å in metal-rich and the O vi 1035Å doublet in metal-poor environments. The C iv emission line were found in some metal-poor starburst galaxies that are good candidates for suppressed or minimal wind signatures (Senchyna et al. 2017;Berg et al. 2019a,b). As discussed by Gray et al. (2019a), these observations might be associated with kinematic features of suppressed bipolar superwinds rather than resonant scattering mentioned by Berg et al. (2019a). An O vi absorption line associated with a weak outflow was found in Haro 11, while the O vi emission luminosity suggests some cooling loss (Grimes et al. 2007). The O vi λ1035 doublet absorption identified in a gravitationally lensed, galaxy has the features of weak low-ionization winds (Chisholm et al. 2018), which may also be explained by the mass-loss effect under non-equilibrium conditions. Moreover, observations of a star-forming galaxy depict an extended halo in the O vi image, and a weak O vi absorbing outflow in the spectrum (Hayes et al. 2016), which might be an indication of suppressed winds.
The enhancements of the C iv and O vi lines could be related to substantial radiative cooling as suggested by Gray et al. (2019a) and Gray et al. (2019b). Time-dependent NEI states calculated by de Avillez & Breitschwerdt (2012) also imply that O vi can be produced in NEI at 10 4.2−5 K below the temperatures that produce O vi in CIE. Similarly, our NEI calculations (before the shell; see Figure 1) also show that the C iv and O vi emission lines do not behave the same in CIE and NEI, especially in the outflow region strongly impacted by radiative cooling.
Time-dependent ionization processes could also be sensitive to time-evolving ionizing sources. Our radiation field was made by Starburst99 for a typical age of 1 Myr. The radiation field calculated by Starburst99 can evolve with the age, which can affect the formation of radiative cooling in starburst-driven superwinds. Our future hydrodynamic simulations with time-evolving radiation fields will help us to understand better the implication of time-dependent non-equilibrium ionization for star-forming regions.