Catastrophic Cooling in Superwinds. II. Exploring the Parameter Space

Galactic-scale ionized gaseous outflows, the so-called "superwinds" as defined by Heckman et al. (1990), are commonly observed in starburst galaxies with extreme star formation activity (Heckman et al. 1987; McCarthy et al. 1987; Fabbiano 1988; Heckman et al. 1990, 1993; Lehnert & Heckman 1995, 1996; Dahlem 1997; Dahlem et al. 1997; Heckman 2002; Rupke et al. 2002, 2005; Martin 2005; Veilleux et al. 2005). It has been understood that superwinds are driven by the kinetic and thermal energies deposited into the interstellar medium (ISM) by mechanical and radiative feedback from young massive stars (Abbott 1982; Leitherer et al. 1992; Hopkins et al. 2012) and supernova (SN) explosions (Thornton et al. 1998; Mac Low & Ferrara 1999; Scannapieco et al. 2002; Creasey et al. 2013), typically in super star clusters (SSCs; Holtzman et al. 1992; O'Connell et al. 1995; Satyapal et al. 1997; Turner et al. 2003; Melo et al. 2005; Smith et al. 2006; Gilbert & Graham 2007; Galliano et al. 2008). These superwinds are found to be very common in intermediate- and high-redshift star-forming galaxies (Pettini et al. 2001, 2002; Wilman et al. 2005; Steidel et al. 2010), so they may play a crucial role in enriching the intergalactic medium (IGM) in the early universe (Nath & Trentham 1997; Ferrara et al. 2000; Lloyd-Davies et al. 2000; Heckman 2002).

The kinetic and thermal energies from OB stars displace the surrounding ISM and heat it up to ∼10⁷ K, which can result in large-scale expanding shells, the so-called "superbubbles" (Castor et al. 1975; Bruhweiler et al. 1980; Mac Low & McCray 1988; Norman & Ikeuchi 1989), and the production of diffuse X-ray emission (Chu & Mac Low 1990; Chu et al. 1995; Magnier et al. 1996; Strickland et al. 2002; Silich et al. 2005; Añorve-Zeferino et al. 2009). Superbubbles have been detected around OB stellar clusters (Cash et al. 1980; Abbott et al. 1981; Mac Low & McCray 1988; Oey & Massey 1995; Dove et al. 2000; Reynolds et al. 2001), as well as in many star-forming galaxies (Veilleux et al. 1994; Heckman et al. 1995; Marlowe et al. 1995; Weiß et al. 1999; Strickland & Stevens 1999; Heckman et al. 2001; Sakamoto et al. 2006; Tsai et al. 2009). These wind-blown bubbles carry off metal-rich material created by massive stars and SN explosions and enrich the ISM and IGM with the recently produced metals. The dynamics of superbubbles have been modeled by several authors (e.g., Castor et al. 1975; Weaver et al. 1977; Koo & McKee 1992a, 1992b; Silich et al. 2005). The theoretical calculations by Weaver et al. (1977) suggested the presence of highly ionized ions such as O vi within the hot interior of the superbubble, where the temperature spans the range ∼ 10⁵–10⁶ K. Moreover, the analytical study by Silich et al. (2005) hinted at the hard (2–8 keV) X-ray emission from the compact hot thermalized ejecta and the soft (0.3–2 keV) diffuse X-ray emission from the superbubble interior.

The mechanical feedback from stellar winds and SNe has been modeled using adiabatic fluid flow by a number of authors (Castor et al. 1975; Weaver et al. 1977; Chevalier & Clegg 1985; Cantó et al. 2000). The adiabatic solutions by Castor et al. (1975) and Weaver et al. (1977) were derived by taking an expanding shell surrounded by the H i and H ii regions of a uniform density that leads to the formation of a bubble. Similarly, the adiabatic steady wind solution by Chevalier & Clegg (1985) was obtained by adopting a continuous freely expanding stationary wind without any ambient medium. The analytical solution of Chevalier & Clegg (1985, hereafter CC85) describes strong superwinds driven by the kinetic and thermal energies supplied by SNe and massive stars of a stellar cluster. The CC85 model provided the first approximate analytic solution to the adiabatic hydrodynamic properties of the freely streaming stationary superwinds, indicating that density, temperature, and thermal pressure of the superwind around the SSC decrease with radius r as ρ_w ∝ r⁻², T_w ∝ r^−4/3, and P_w ∝ r^−10/3, respectively, while the wind velocity u_w approaches the adiabatic terminal speed V_A,∞. The numerical simulations of the CC85 model conducted by Cantó et al. (2000) demonstrated that both the analytic adiabatic solutions and their numerical calculations are fully consistent for spherically symmetric homogeneous winds.

Moreover, the semianalytic numerical studies by Silich et al. (2003) and Silich et al. (2004) explored the impact of strongly radiative cooling, called "catastrophic cooling," on density, temperature, and velocity of an outflow. In particular, Silich et al. (2004) found that the wind temperature predicted by the radiative solution has a departure from the adiabatic approximation of T_w ∝ r^−4/3 in massive compact clusters. Moreover, an increase in the ejected gas metallicity Z was found to enhance radiative cooling within the cluster volume (Silich et al. 2004; Tenorio-Tagle et al. 2005). The semianalytic modeling and simulations by Tenorio-Tagle et al. (2007) estimated the locus for catastrophic cooling in the plane of cluster size (R_sc) versus the mechanical luminosity (${L}_{\mathrm{mech}}=\tfrac{1}{2}{\dot{M}}_{\mathrm{sc}}{V}_{A,\infty }^{2}$), which is a function of the total stellar mass-loss rate (${\dot{M}}_{\mathrm{sc}}$) and the adiabatic terminal speed (V_A,∞).

Recent observations are reported to exhibit signs of catastrophic cooling and the suppression of superwinds and superbubbles in SSCs in several starburst galaxies, such as M82 (Smith et al. 2006; Westmoquette et al. 2014), NGC 2366 (Oey et al. 2017), and NGC 5253 (Turner et al. 2017), where the SSCs are embedded within 1–2 pc, dense, ultracompact, high-pressure gas. Similarly, the most extreme Green Peas (GPs; Jaskot et al. 2017) show kinematics consistent with suppressed superwinds, contrary to other starbursts. Recent ALMA observations of distant star-forming galaxies reveal the presence of extended [C ii] halos that could also be produced by starburst-driven cooling outflows (Fujimoto et al. 2019, 2020; Pizzati et al. 2020). These starburst-driven strongly cooled outflows could be explained by the radiative solution semianalytically derived by Silich et al. (2003, 2004). More recently, Gray et al. (2019a) investigated the occurrence of catastrophic cooling through hydrodynamic simulations including radiative nonequilibrium cooling, which showed the strong enhancement of highly ionized ions such as C iv and O vi.

In the recent work by Gray et al. (2019a, Paper I), the hydrodynamic simulations with the nonequilibrium atomic chemistry and cooling package MAIHEM (Gray et al. 2015; Gray & Scannapieco 2016; Gray et al. 2019b) were performed assuming a default solar value for the metallicity Z of the outflow and a typical mass-loss rate of 10⁻² M_⊙ yr⁻¹. As found by Silich et al. (2004), Tenorio-Tagle et al. (2005), and Gray et al. (2019b), the metallicity can have a significant role in strongly enhancing radiative cooling. Moreover, the mechanical luminosity, i.e., the mass-loss rate and wind terminal speed, can affect the domain of catastrophic cooling (Tenorio-Tagle et al. 2007).

In this paper, we therefore investigate a larger parameter space to explore the effect of metallicity Z, mass-loss rate ${\dot{M}}_{\mathrm{sc}}$, and wind terminal velocity V_∞ on the behavior of catastrophic cooling and superwind suppression. Although our models are parameterized in terms of ${\dot{M}}_{\mathrm{sc}}$ and V_∞, we note that these are the effective values that may result after mass loading or heating efficiency effects. As in Paper I, we use one-dimensional, spherically symmetric hydrodynamic simulations coupled to radiative thermal functions using the nonequilibrium atomic chemistry and cooling package MAIHEM (Gray et al. 2019a) to run an extensive grid of models. In addition, we also investigate some emission-line diagnostics for these conditions. These are generated by post-processing the MAIHEM models with the CLOUDY photoionization code.

This paper is organized as follows. Section 2 describes the configuration of our hydrodynamic simulations of superwinds, including the radiative thermal functions and initial and boundary conditions. Section 3 describes the implementation of the photoionization modeling, as well as computing population synthesis models for SSCs. The hydrodynamic results from our MAIHEM simulations are presented in Section 4, followed by the emission lines calculated by CLOUDY and diagnostic diagrams in Section 5. Our results are compared with observations in Section 6 and are discussed in Section 7. Our conclusions are given in Section 8.

We solve the hydrodynamic equations coupled to the radiative cooling and heating functions using the nonequilibrium atomic chemistry and cooling package MAIHEM (Models of Agitated and Illuminated Hindering and Emitting Media; Gray et al. 2015; Gray & Scannapieco 2016; Gray et al. 2019b), which is a modified version of the adaptive mesh hydrodynamics code FLASH v4.5 (Fryxell et al. 2000). We use the directionally unsplit pure hydrodynamic solver (Lee et al. 2009; Lee & Deane 2009; Lee 2013), together with the second-order Monotone Upwind-centered Scheme for Conservation Laws (MUSCL)–Hancock scheme (van Leer 1979) developed for arbitrary Lagrangian Eulerian (ALE) and smoothed particle hydrodynamics (SPH) methods. ² To avoid odd–even instabilities, we also employ a hybrid type of the Riemann solver (Toro et al. 1994), combining the Roe solver (Roe 1981) and the Harten–Lax–van Leer–Einfeldt (HLLE) solver (Einfeldt 1988; Einfeldt et al. 1991). The HLLE solver is a modified version of the HLL solver that is positively conservative under certain wave-speed conditions in order to provide more stable and diffusive solutions in strongly shocked regions.

Following Silich et al. (2004), the equations of hydrodynamics applicable to spherically symmetric superwinds in a steady state with the radiative cooling for a spherically symmetric SSC of radius R_sc without accounting for the gravitational attraction from the SSC are

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{d}{{dr}}\left({\rho }_{w}{u}_{w}{r}^{2}\right)={q}_{m},\end{eqnarray} \tag{ 1 }$$

$$\begin{eqnarray}&&{\rho }_{w}{u}_{w}\displaystyle \frac{{{du}}_{w}}{{dr}}+\displaystyle \frac{{{dP}}_{w}}{{dr}}=-{q}_{m}{u}_{w},\end{eqnarray} \tag{ 2 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{d}{{dr}}\left[{\rho }_{w}{u}_{w}{r}^{2}\left(\displaystyle \frac{{u}_{w}^{2}}{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{{P}_{w}}{{\rho }_{w}}\right)\right]={q}_{e}-{q}_{c}+{q}_{h},\end{eqnarray} \tag{ 3 }$$

where r is the radial distance, u_w the wind velocity, ρ_w the wind density, P_w the wind thermal pressure, ${q}_{m}=3{\dot{M}}_{\mathrm{sc}}/4\pi {R}_{\mathrm{sc}}^{3}$ the total mass-loss rate per unit volume, ${q}_{e}=3{\dot{E}}_{\mathrm{sc}}/4\pi {R}_{\mathrm{sc}}^{3}$ the total energy deposition rate per unit volume, ${\dot{M}}_{\mathrm{sc}}$ the total mass-loss rate, ${\dot{E}}_{\mathrm{sc}}$ the total energy deposition rate, γ = 5/3 the ratio of specific heats, ${q}_{c}={n}_{w}^{2}{\rm{\Lambda }}(Z,{T}_{w})$ the cooling rate per unit volume, q_h = n_w Γ(Z) the heating rate per unit volume, n_w the wind number density, Λ(Z, T_w ) the cooling function (Raymond et al. 1976; Oppenheimer & Schaye 2013), Γ(Z) the heating function (see, e.g., Wolfire et al. 2003; Gnedin & Hollon 2012; Gray & Scannapieco 2016; Gray et al. 2019b), T_w the wind temperature, and Z the metallicity. We adopt the gamma-law equation of state for an ideal gas, P_w = (γ − 1)ρ_{w w} , where _w is the internal energy. This accounts for changes in the mean atomic mass due to evolution of the ionization states of all the species. Outside the SSC (r > R_sc), q_m and q_e vanish (CC85). The wind model assumes that radiative thermal effects within the SSC are negligible, so q_c and q_h vanish at r < R_sc. Taking q_c = 0 and q_h = 0 at r > R_sc, we recover the adiabatic solution in CC85 and Cantó et al. (2000). The wind terminal velocity is also defined as ${V}_{\infty }={(2{\dot{E}}_{\mathrm{sc}}/{\dot{M}}_{\mathrm{sc}})}^{1/2}$.

Assuming the wind of a uniform density, solving Equations (1)–(3) yields the wind density and pressure at r > R_sc away from the SSC (Silich et al. 2004):

$$\begin{eqnarray}&&{\rho }_{w}=\displaystyle \frac{{\dot{M}}_{\mathrm{sc}}}{4\pi {u}_{w}{r}^{2}},\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}&&\displaystyle \frac{{{dP}}_{w}}{{dr}}=\displaystyle \frac{{\dot{M}}_{\mathrm{sc}}}{4\pi {r}^{3}}\displaystyle \frac{(\gamma -1)r({q}_{h}-{q}_{c})-2\gamma {u}_{w}{P}_{w}}{{\rho }_{w}({u}_{w}^{2}-{c}_{s}^{2})},\end{eqnarray} \tag{ 5 }$$

where ${c}_{s}={(\gamma {P}_{w}/{\rho }_{w})}^{1/2}={(\gamma {k}_{{\rm{B}}}{T}_{w}/\mu )}^{1/2}$ is the sound speed, μ the mean mass per wind particle, and k_B the Boltzmann constant. Following the stationary wind solution of Silich et al. (2004), at r = R_sc, the Mach number M = u_w /c_s = 1 (CC85), so the wind velocity u_w is equal to the sound speed c_s . Following CC85, the wind velocity is u_w = V_∞/2 at the cluster boundary r = R_sc. Equation (4) implies that the wind density at the cluster boundary should be ${\rho }_{w}\,=\,{\dot{M}}_{\mathrm{sc}}/2\pi {V}_{\infty }{R}_{\mathrm{sc}}^{2}$. As c_s = u_w at r = R_sc, the wind temperature and pressure at the cluster boundary should be ${T}_{w}=\mu {V}_{\infty }^{2}/4\gamma {k}_{{\rm{B}}}$ and ${P}_{w}={\rho }_{w}{V}_{\infty }^{2}/4\gamma$, respectively. These values of the wind density and temperature are used as the boundary conditions at r = R_sc in our hydrodynamic simulations as described in the following subsection.

2.1. Initial and Boundary Conditions

The hydrodynamic simulations are computed from t = 0 until an age of 1 Myr at checkpoint intervals of 0.1 Myr. This final age is long enough to generate a mature outflow but too short for any appreciable evolution in the spectral energy distribution (SED) and stellar population. It describes typical superwinds driven by extremely young starbursts. The hydrodynamic model employs one-dimensional, spherical geometry with the initial radius equal to the SSC R_sc, which is set to 1 pc in our study. This corresponds to the flow injection radius. The outer boundary is set to 250 pc in our hydrodynamic simulations, which covers the expansion of fast winds with high mass-loss rates at their final age. We note that the associated Strömgren spheres are often larger, and the outer boundaries in our CLOUDY models are unlimited in radius (described in Section 3). The hydrodynamic models are simulated on a base grid consisting of N_x = 512 blocks that are allowed up to two levels of adaptive mesh refinement. This gives a maximum resolution of 0.244 pc. The adaptive mesh algorithm (paramesh; MacNeice et al. 2000) employed by FLASH uniformly covers the physical space, whose refinement criterion is an error estimator (Löhner 1987) that is configured to follow the changes in density, temperature, pressure, and velocity. The ambient medium surrounding the SSC does not have any initial kinematics (u_w = 0).

The initial conditions of the wind velocity, temperature, and density at the cluster boundary for the hydrodynamic simulations are adopted based on the analytic radiative solutions (Silich et al. 2004), which are similar to the CC85 adiabatic solutions. The initial values of the temperature and density in the surrounding regions are set to the temperature T_amb—which is determined by photoionization of the SSC—and density n_amb of the ambient medium, respectively. The initial wind density ρ_w,sc and temperature T_w,sc at the cluster boundary (r = R_sc) are then calculated by assuming the analytic solution:

$$\begin{eqnarray}&&{\rho }_{w,\mathrm{sc}}=\displaystyle \frac{{\dot{M}}_{\mathrm{sc}}}{4\pi {R}_{\mathrm{sc}}^{2}{u}_{w,\mathrm{sc}}},\end{eqnarray} \tag{ 6 }$$

$$\begin{eqnarray}&&{T}_{w,\mathrm{sc}}=\displaystyle \frac{1}{\gamma }\displaystyle \frac{\mu }{{k}_{{\rm{B}}}}{u}_{w,\mathrm{sc}}^{2},\end{eqnarray} \tag{ 7 }$$

where u_w,sc is the wind velocity at r = R_sc given by the wind terminal speed V_∞ (Chevalier & Clegg 1985; Cantó et al. 2000),

$$\begin{eqnarray}&&{u}_{w,\mathrm{sc}}=\tfrac{1}{2}{V}_{\infty }.\end{eqnarray} \tag{ 8 }$$

We initialize the boundary wind velocity to V_∞/2. This initial wind velocity configuration is based on the wind terminal speed that easily allows us to compare the radiative solution to the adiabatic solution. For the effective area of the outflow ${\rm{\Omega }}{R}_{\mathrm{sc}}^{2}$ in Equation (6), we use Ω = 4π corresponding to a fully isotropic outflow, while Gray et al. (2019a) adopted Ω = π describing an anisotropic broken outflow extended perpendicularly to the galactic plane.

Following the CC85 adiabatic solutions, the wind velocity at r = R_sc is assumed to be equal to the local sound speed, u_w,sc = c_s,sc, which corresponds to a Mach number of 1. The wind terminal speed V_∞ is related to the mechanical luminosity as

$$\begin{eqnarray}&&{L}_{\mathrm{mech}}=\tfrac{1}{2}{\dot{M}}_{\mathrm{sc}}{V}_{\infty }^{2}.\end{eqnarray} \tag{ 9 }$$

After setting the initial temperature and density from the assumed ambient density, the thermal pressure and internal energy are determined by solving Equations (1)–(3) in MAIHEM, using the equation of state. Thus, the wind solution develops from these initial conditions.

The initial thermal pressure derived from the hydrodynamic simulation is found to be consistent with the analytically derived initial pressure:

$$\begin{eqnarray}&&{P}_{w,\mathrm{sc}}=\displaystyle \frac{1}{\gamma }\displaystyle \frac{{\dot{M}}_{\mathrm{sc}}{u}_{w,\mathrm{sc}}}{4\pi {R}_{\mathrm{sc}}^{2}}.\end{eqnarray} \tag{ 10 }$$

Alternatively, we could employ the initial analytic pressure P_w,sc and density ρ_w,sc and determine the initial temperature T_w,sc by solving the equation of state in MAIHEM that would provide the same analytic temperature given in Equation (7). The internal energy per unit mass can also be calculated from the internal pressure and density as follows:

$$\begin{eqnarray}&&{\epsilon }_{w,\mathrm{sc}}=\displaystyle \frac{1}{(\gamma -1)}\displaystyle \frac{{P}_{w,\mathrm{sc}}}{{\rho }_{w,\mathrm{sc}}},\end{eqnarray} \tag{ 11 }$$

which is calculated by solving the equation of state in MAIHEM. The internal and kinetic energies per unit mass are used to initialize the total energy per unit mass at the cluster boundary (${E}_{w,\mathrm{sc}}={\epsilon }_{w,\mathrm{sc}}+\tfrac{1}{2}| {u}_{w,\mathrm{sc}}{| }^{2}$).

Following Gray et al. (2019a, 2019b), the gas ejected by the cluster at R_sc is approximated to be in collisional ionization equilibrium (CIE) at the boundary, so it depends only on the initial wind temperature and ionizing UV background. The photoionization code CLOUDY is called to generate the CIE ionization fractions that are a function of temperature and applicable ionizing UV background. Our procedure makes tabulated grids of ionization fractions for all the ions in CIE that span the temperature range of 10²–10⁸ K, which are linearly interpolated over to yield the ionization fractions at the outflow boundary for a given temperature during our hydrodynamic simulation. The tabulated grids incorporate all the processes considered by CLOUDY, including collisional excitations, recombination, and collisional ionization. The ionization fractions of the ejected gas at the outflow boundary are then calculated by MAIHEM at runtime through an interpolation on the CIE tabulated grids generated by CLOUDY.

Table 1 summarizes the parameters used in our MAIHEM hydrodynamic simulations, including the metallicity used in MAIHEM simulations (Z), the mass-loss rate (${\dot{M}}_{\mathrm{sc}}$), and the wind terminal speed (V_∞) in the first three columns. It also lists the metallicities used in Starburst99 models for the Geneva-Rot stellar evolution, high-resolution spectra, and UV spectrum (Columns (4)–(6)). The last three columns present Starburst99 outputs: the total luminosity, L_tot, the fraction of ionizing photons in the H i continuum relative to the total luminosity, f(H⁺), and the ionizing luminosity, L_ion. We consider the gas metallicity of $\hat{Z}\equiv Z/$ Z_⊙ = 1, 0.5, 0.25, and 0.125, where Z/Z_⊙ = 1 corresponds to the ISM abundances listed in Table 2. We adopt the typical ISM abundances of Savage et al. (1977), together with the ISM gas-phase oxygen abundance of Meyer et al. (1998), as the baseline solar metallicity (see Table 2). Typical ISM grains with a dust-to-metal mass ratio of M_d/M_Z = 0.2 are incorporated into our initial models of the ionization fraction computed by CLOUDY (see also Section 3). Dust grains are not included as a distinct species in the MAIHEM simulations, but we use C/O ratios that account for depletion as described in Section 3. For the fiducial models with the solar metallicity (Z/Z_⊙ = 1), we calculate the boundary density ρ_w,sc and temperature T_w,sc using the total mass-loss rate ${\dot{M}}_{\mathrm{sc}}$ of 10⁻¹, 10⁻², 10⁻³, and 10⁻⁴ M_⊙ yr⁻¹ and the wind terminal speed V_∞ of 250, 500, and 1000 km s⁻¹. For the models with subsolar metallicities, namely, Z/Z_⊙ = 0.5 , 0.25, and 0.125, we scale the mass-loss rates and wind velocities of the solar models according to ${\dot{M}}_{\mathrm{sc}}\propto {Z}^{0.72}$ (Mokiem et al. 2007) and V_∞ ∝ Z^0.13 (Vink et al. 2001), since we are interested only in extremely young, pre-SN SSCs that are dominated by stellar winds.

Table 1. Summary of Parameters in Hydrodynamic Simulations by MAIHEM and Stellar Population Models by Starburst99

Hydro	MAIHEM Outflow Parameters		Starburst99 Inputs			Starburst99 Output
Z	${\dot{M}}_{\mathrm{sc}}(\propto {Z}^{0.72})$	V∞( ∝ Z0.13)	Z	Z	Z	$\mathrm{log}{L}_{\mathrm{tot}}$	f(H+)	$\mathrm{log}{L}_{\mathrm{ion}}$
(Z⊙)	(M⊙ yr−1)	(km s−1)	Geneva-Rot	High Res. Spec.	UV Spec.	(erg s−1)		(erg s−1)
1.0	1.0 × [10−1, 10−2, 10−3, 10−4]	[250, 500, 1000]	0.014	0.020	Solar	42.97	0.41	42.58
0.5	0.607 × [10−1, 10−2, 10−3, 10−4]	[229, 457, 914]	0.008	0.008	LMC/SMC	42.97	0.42	42.60
0.25	0.369 × [10−1, 10−2, 10−3, 10−4]	[209, 418, 835]	0.008	0.008	LMC/SMC	42.97	0.42	42.60
0.125	0.224 × [10−1, 10−2, 10−3, 10−4]	[191, 382, 736]	0.002	0.008	LMC/SMC	42.98	0.50	42.67

Note. The ambient temperature (T_amb) is set to the value calculated by the cloudy model for the given Starburst99 SED. Other MAIHEM parameters are as follows: the cluster radius R_sc = 1 pc, the ambient density n_amb = 1, 10, ..., 10³ cm⁻³, the number of blocks N_x = 512, and time intervals t = 0.1, 0.2, ..., 1.0 Myr. Other Starburst99 parameters are as follows: the total stellar mass M_⋆ = 2.05 × 10⁶ M_⊙, and age t = 1 Myr.

Download table as:  ASCIITypeset image

Table 2. ISM Abundances and Depletion Factors for CLOUDY Models

Element	n(X)/n(H)	(1 − fdpl(X))
He	0.098	1.0
C	2.87 × 10−4	0.4
N	7.94 × 10−5	1.0
O	3.19 × 10−4	0.6
Ne	1.23 × 10−4	1.0
Na	3.16 × 10−7	0.2
Mg	1.26 × 10−5	0.2
Si	3.16 × 10−6	0.03
S	3.24 × 10−5	1.0
Ar	2.82 × 10−6	1.0
Ca	4.10 × 10−10	10−4
Fe	6.31 × 10−7	10−2

Note. The abundance n(X)/n(H) of the element X is by number relative to the hydrogen abundance. The depletion factor f_dpl(X) specifies the depletion fraction of the element X onto dust grains. The carbon abundance is set in a way to provide [n(C)(1 − f_dpl(C)]/[n(O)(1 − f_dpl(O)] = 0.6, described in text.

Download table as:  ASCIITypeset image

For the ambient medium in our MAIHEM simulations, we adopt the photoionized temperature structure created by the ionizing SED from the parent SSC. This is done using a single CLOUDY model applied to the density profile predicted by MAIHEM from a preliminary simulation, for the fiducial age of 1 Myr, described in detail in Section 3.1 (the "PI model"). This sets the ionized radius for the final MAIHEM model and determines how long the outflow expands into an ionized medium. For most models, the ambient medium remains ionized at 1 Myr. For optically thick models transitioning to a neutral ISM (see Section 5), which are mostly those with n_amb = 10³ cm⁻³, we adopt an ambient temperature of T_amb = 5 × 10³ K. The ionizing SED is generated by the Starburst99 population synthesis model described in Section 3 for the fiducial model of M_⋆ = 2.05 × 10⁶ M_⊙ at age 1 Myr. This SED is also included as a background UV radiation field in the MAIHEM simulation by using a simple r⁻² decrease in flux. MAIHEM does not calculate detailed radiative transfer, but since the wind is low density and therefore optically thin, this is a reasonable assumption. Initially, the entire surrounding domain outside the SSC is assumed to have ambient densities of n_amb = 1, 10, 10², and 10³ cm⁻³.

2.2. Nonequilibrium Cooling Functions

To include hydrodynamic heating and background UV radiation, the cooling and heating functions of atomic species are calculated using nonequilibrium ionization (NEI) conditions. The cooling rate per volume q_c (Z, T_w ) is determined by the MAIHEM cooling routine using the cooling function Λ(Z, T_w ) based on the gas metallicity Z, in addition to the gas temperature T_w found by hydrodynamic solutions for given physical conditions and background radiation. We use the cooling routine implemented by Gray et al. (2015), which extended the ion-by-ion cooling efficiencies Λ_i (T_w ) of Gnat & Ferland (2012) down to 5000 K, for the given ions. Additionally, Gray et al. (2015) included NEI in MAIHEM, which can track NEI and recombination of species. Furthermore, as fully described by Gray et al. (2019b), the MAIHEM package takes account of the heating rate per volume, q_h (Z), by calculating the heating function Γ(Z) for the outflow of metallicity Z. This atomic chemistry was implemented by incorporating collisional ionization rates (Voronov 1997), radiative recombination rates (Badnell 2006), and dielectronic recombination rates (see Table 1 in Gray et al. 2015). With the inclusion of the atomic chemistry package, our simulations are performed using the multispecies extension for the equation of state that accounts for the change in average atomic weight on the thermodynamic quantities. The species and nonequilibrium cooling were incorporated into MAIHEM. These two capabilities, namely, NEI and radiative thermal function, allow us to use MAIHEM to investigate the radiative solution of the superwind model.

The package MAIHEM computes the photoionization and photoheating rate for each ion (Verner & Yakovlev 1995; Verner et al. 1996) to include the effect of an ionizing UV background (Gray & Scannapieco 2016). The most recent version of MAIHEM updated by Gray et al. (2019b) is able to track NEI and recombination of 84 species across 13 elements, namely, hydrogen (H i–H ii), helium (He i–He iii), carbon (C i–C vii), nitrogen (N i–N viii), oxygen (O i–O ix), neon (Ne i–Ne xi), sodium (Na i–Na vi), magnesium (Mg i–Mg vi), silicon (Si i–Si vi), sulfur (S i–S vi), argon (Ar i–Ar vi), calcium (Ca i–Ca vi), and iron (Fe i–Fe vi). The cooling efficiencies of the new ions in the expanded network of Gray et al. (2019b) were computed using CLOUDY, in addition to the method used by Gray et al. (2015). The tabulated grids of cooling efficiencies are read into MAIHEM during the initialization for each model. The latest-version models include the column densities of N v and O vi (Gray et al. 2019b).

As mentioned in Gray et al. (2019a), the cooling efficiency Λ_i (T_w ) of species i, the number density n_i of species i, and the number of electrons n_e yield the cooling rate per volume:

$$\begin{eqnarray}&&{q}_{c}(Z,{T}_{w})={n}_{w}^{2}{\rm{\Lambda }}(Z,{T}_{w})=\sum _{i}{n}_{i}{n}_{e}{{\rm{\Lambda }}}_{i}({T}_{w}),\end{eqnarray} \tag{ 12 }$$

which is used in the computations of fluid models to provide radiative cooling. As the ion-by-ion cooling efficiencies Λ_i (T_w ) are calculated based on the temperature determined by the hydrodynamic simulations for specified physical conditions and background radiation, we see how the mechanical feedback (via ${\dot{M}}_{\mathrm{sc}}$ and V_∞) can lead to the catastrophic cooling regime. Furthermore, the cooling function Λ(Z, T_w ) is a function of the metallicity specified for the ejected gas.

The heating rate per volume is calculated using the photoheating efficiency Γ_i and number density n_i of species i as follows (Gray et al. 2019b):

$$\begin{eqnarray}&&{q}_{h}(Z)={n}_{w}{\rm{\Gamma }}(Z)=\sum _{i}{n}_{i}{{\rm{\Gamma }}}_{i},\end{eqnarray} \tag{ 13 }$$

where the ion-by-ion photoheating efficiencies Γ_i are estimated by MAIHEM (see Gray & Scannapieco 2016) using a background UV SED and the photoionization cross section of species i, taken from Verner & Yakovlev (1995) and Verner et al. (1996).

We next carry out photoionization modeling to calculate emissivities of a set of UV/optical emission lines that are typically used for diagnostic purposes. These are used to build diagnostic diagrams (Section 5.1), which we compare with observations (Section 6). The photoionization modeling is implemented as post-processing using the program CLOUDY v17.02, which is described by Ferland et al. (1998, 2013, 2017). The outflow temperature T_w and density n_w as a function of the radius r obtained from our hydrodynamic simulations are used as inputs to our photoionization models. Additionally, the SED of the SSC generated by Starburst99 (Leitherer et al. 1999, 2014), together with its predicted ionizing luminosity, as described in Section 3.1, is used as an ionizing source for our photoionization models. The stellar evolution and atmosphere models employed by Starburst99 use different increments in metallicity, and so the metallicity in each of our photoionization models is set to the closest match with the metallicity Z (i.e., Z/Z_⊙ = 1, 0.5, 0.25, and 0.125) used in the corresponding hydrodynamic simulation, as shown in Table 1.

Similar to Gray et al. (2019a), we model the ionization states for two cases using the same hydrodynamically generated density distribution from MAIHEM: a purely photoionized model based on resetting the entire gas distribution to neutral temperature, and therefore without any hydrodynamic thermal contributions (PI case); and a model with both photoionization and hydrodynamic collisional ionization (CPI case), based on the MAIHEM temperature distribution. These are analogous to models B and C in Paper I. We use a single CLOUDY model to generate a more realistic ionization structure, although it is necessarily CIE, since CLOUDY cannot model non-CIE. We summarize these models as follows (Figure 1):

• 1.  
  PI case: This model is only used to generate inputs for the CPI case. The ambient temperature in MAIHEM is initially set to an isothermal value of 1000 K to estimate initial density profiles. Together with the ionizing SED and luminosity from Starburst99, they are employed by CLOUDY to calculate the photoionization and emissivities due to the stellar radiation of the SSC. The CLOUDY calculations stop when the temperature drops below 950 K. The resulting temperature profile in the PI case is then used to obtain the photoionized radius and mean ambient temperature for the ionized portion of the ambient medium and shell in the CPI case (see below and Section 2.1).
• 2.  
  CPI case: Both the density and temperature profiles generated by MAIHEM, together with the ionizing SED and luminosity produced by Starburst99, are used by CLOUDY to compute the ionization states and emissivities resulting from both the stellar radiation of the SSC and thermal contributions from the hydrodynamic outflow (see Figure 1). The temperature profile found in the PI case is applied to the ambient medium and defines the radius at which the model transitions to neutral and becomes optically thick. The mean ambient ionized temperature in the PI case is adopted for running the MAIHEM model. It is also adopted for the temperature in the ambient medium surrounding the shell, as well as the ionized, isothermal region of the shell, roughly ∼1–2 pc outward from the shell inner boundary. The CPI calculations stop according to the optical depth (H⁺) defined by the PI model.

Zoom In Zoom Out Reset image size

Thus, the PI CLOUDY model is used to define the ambient temperature profile in the ambient medium for the hydrodynamic simulation by MAIHEM (see Section 2), as well as the ambient medium in the CPI CLOUDY model. The MAIHEM simulation produces the NEI states and temperature profile of the outflow. However, we only employ the temperature profile generated by MAIHEM for the CPI model, and not the ionization states, so our CLOUDY models calculate only CIE conditions. NEI conditions, which can be estimated by including both the ionization states and temperature structure predicted by MAIHEM, generate higher fluxes in certain highly ionized lines such as C iv and O vi, as shown in Paper I. Our CPI models, which combine the thermal effects and density distribution from both the hydrodynamic outflow and photoionization, are the final, realistic models.

For the CLOUDY models with Z/Z_⊙ = 1, we adopt the ISM abundances derived by Savage et al. (1977) and the ISM oxygen abundance by Meyer et al. (1998), together with their metal depletion onto dust grains (see Table 2), the typical ISM grains, and the C/O ratio parameterized by the metallicity. Several observations implied that the C/O ratio depends on the metallicity (e.g., Garnett et al. 1995, 1999). Considering the C/O measurements by Garnett et al. (1999), the C/O ratio correlates with the O/H abundance as logC/O = 3.84 + 1.16logO/H, so we get C/O = 0.16 for the Small Magellanic Cloud (SMC; logO/H = −4, i.e., Z/Z_⊙ ∼ 1/5), C/O = 0.46 for the Large Magellanic Cloud (LMC; logO/H = −3.6, i.e., Z/Z_⊙ ∼ 1/2), and C/O = 0.6 for the local Galactic ISM (logO/H = −3.5, i.e., Z/Z_⊙ ∼ 1). The carbon abundance is configured to yield the depleted C/O ratio associated with the relevant metallicities (see Table 2). We apply the CLOUDY "metals" command to the models with the metallicities Z/Z_⊙ = 0.5, 0.25, and 0.125 to scale down the baseline solar abundances of elements heavier than helium. The same metal depletion and adjusted C/O ratios are also incorporated into CLOUDY models used by MAIHEM for initializing CIE states of the ejected gas, as well as for metal species used in the NEI calculations by our MAIHEM simulations.

For our CLOUDY models, we also include the typical ISM grains with a dust-to-metal mass ratio of M_d/M_Z = 0.2, which is typically associated with evolved galaxies (De Vis et al. 2019). Previously, the dust-to-gas versus metallicity correlation (see e.g., Hirashita 1999; Edmunds 2001) corresponds to M_d/M_Z ∼ 0.5. However, the recent study by De Vis et al. (2019) suggested that evolved galaxies with metallicities above logO/H = −3.8 have a more or less constant dust-to-metal ratio of M_d/M_Z ≈ 0.214, while significantly lower values are expected for unevolved galaxies (logO/H < −3.8). As we do not consider implications of dust grains in this paper, we assume the standard ISM grains having a constant ratio of M_d/M_Z = 0.2. The dust grains, along with depleted metallic species, are included in CLOUDY models for initializing CIE states of the outflow in MAIHEM, but not in the NEI module.

3.1. Population Synthesis Models

In this section, we present the results of our hydrodynamic simulations performed by MAIHEM using the parameters listed in Table 1.

To illustrate different features of the temperature and density profiles in the adiabatic and radiative cooling modes, in Figure 2 we plot the temperature T_w and density n_w as a function of radius r. Considering the hydrodynamic simulations displayed as an example figure at t = 1 Myr with the metallicity Z/Z_⊙ = 0.5, mass-loss rate ${\dot{M}}_{\mathrm{sc}}=6.07\times {10}^{-4}$ M_⊙ yr⁻¹, cluster radius R_sc = 1 pc, ambient density n_amb = 10 cm⁻³, and ambient temperature T_amb ∼ 1.4 × 10⁴ K determined by CLOUDY, there are two with negligible radiative cooling (V_∞ = 457 and 914 km s⁻¹) and one with considerable radiative cooling (V_∞ = 229 km s⁻¹). The hydrodynamic models (V_∞ = 457 and 914 km s⁻¹) without radiative cooling are not continuous, freely expanding CC85 winds owing to the presence of a shocked-wind region described by Weaver et al. (1977), but they possess the CC85 adiabatic solution before the bubble in their temperature and density profiles. The hydrodynamic model (V_∞ = 229 km s⁻¹) with strongly radiative cooling does not follow the CC85 adiabatic cooling flow before the shell, and it has a radiative solution more similar to that described by Silich et al. (2004).

The adiabatic cooling mode in V_∞ = 914 km s⁻¹ shown in the example image of Figure 2 has four different regions described in Weaver et al. (1977): (a) freely expanding supersonic wind, (b) shocked-wind region, (c) a narrow shell of shocked ambient medium, and (d) the ambient medium. We see that the temperature and density profiles in region (a) follow the prediction by the CC85 analytic adiabatic solution (dashed line). As the wind propagates into the ambient medium, a hot bubble of shocked wind results, with sharply increased temperature. This leads to the formation of a dense shell due to the confrontation with the ambient medium. The shell has a density higher than the ambient density.

The nonadiabatic, radiative cooling mode is shown in an example figure of Figure 3: this model (V_∞ = 250 km s⁻¹) has all the regions described in the adiabatic cooling mode, except for the shocked-wind region (b) corresponding to the hot bubble. The animated version of this figure shows the evolution for t = 0.1, 0.2,...,1 Myr (available in the online journal). ³ The strong cooling effects suppress the formation of a powerful shocked-wind region and an intense temperature bubble. Moreover, the temperature of the freely expanding supersonic wind does not follow the analytic adiabatic solution (dashed line). It can be seen that radiative cooling strongly reduces the temperature of the freely expanding supersonic wind that leads to the deposition of much lower kinetic and thermal energy in the surrounding ambient medium. Thus, we identify three regions in the strongly radiative cooling (catastrophic cooling) mode: freely expanding supersonic wind, a broad shell with density slightly higher than the ambient density, and the surrounding medium. This model corresponds to a momentum-driven shell.

Figure 3 depicts adiabatic cooling with V_∞ = 1000 km s⁻¹ that does not develop any bubble at ages t ≤ 0.9 Myr (see its animation in Supplementary Material). However, this model generates a bubble at t = 1 Myr, so the hot bubble characteristic of the adiabatic solution can take time to develop. In our grid, the model for V_∞ = 500 km s⁻¹, n_amb = 1000 cm⁻³, Z/Z_⊙ = 1 at ${\dot{M}}_{\mathrm{sc}}={10}^{-3}$ M_⊙ yr⁻¹ generates bubbles at 1.1 Myr, which is greater than the 1 Myr age adopted for all the models. It is important to note that the lack of a hot bubble in young objects does not necessarily indicate that the system is strongly cooling or otherwise not adiabatic. Conversely, the hot, shocked-wind region can also be present in our simulations at t = 1 Myr with radiative cooling (e.g., V_∞ = 229 km s⁻¹ in Figure 2), whose freely expanding wind experiences mildly cooling, but not strong enough to suppress the shocked wind and hot bubble.

The temperature and density profiles of the superwinds simulated by our hydrodynamic models are therefore classified according to the presence of adiabatic/radiative cooling and the bubble as follows:

• 1.  
  Adiabatic wind (AW) mode: The temperature profile is roughly of freely expanding supersonic winds (CC85) that closely follow the adiabatic solution, and it has not yet produced any bubble, which may develop at a later time. The temperature T_w of the expanding wind (region (a) in Figure 2) produced by MAIHEM is within 0.75 and 1.25 times the temperature T_adi analytically derived from the adiabatic solution.
• 2.  
  Adiabatic bubble (AB) mode: This is the classic bubble model (Weaver et al. 1977), where the temperature profile is roughly of adiabatic winds (slightly or without radiative cooling) surrounded by a shocked-wind region (hot bubble; with a thickness >0.73 pc and a mean temperature >10⁵ K). The temperature profile closely follows the adiabatic solution of freely expanding supersonic winds (CC85). The temperature T_w of the expanding wind computed by MAIHEM is within 0.75 and 1.25 of the analytic adiabatic temperature T_adi over the expanding wind region. The hot bubble (region (b) in Figure 2) drives the shell expansion.
• 3.  
  Adiabatic, pressure-confined (AP) mode: The expanding wind region is adiabatic, but the bubble expansion is stalled by high ambient pressure from the dense environment. This corresponds to the pressure-confined model predicted by Silich et al. (2007; see also Koo & McKee 1992a). The temperature profile is roughly of freely expanding supersonic winds (CC85). The expanding wind temperature T_w computed by MAIHEM is between 0.75 and 1.25 times the temperature T_adi analytically derived from the adiabatic solution.
• 4.  
  Catastrophic cooling (CC) mode: The temperature profile is of freely expanding winds with strongly radiative cooling (Silich et al. 2003, 2004) without any noticeable shock regions and weak/no bubble. The expanding wind temperature T_w of the freely expanding wind (region (a)) simulated by MAIHEM is at least 25% lower than the temperature T_adi analytically derived from the adiabatic solution.
• 5.  
  Catastrophic cooling bubble (CB) mode: The temperature profile has strongly radiative cooling but still has the bubble and shocked-wind regions. The expanding wind temperature T_w of the region (a) is at least 25% below the adiabatic temperature T_adi from the analytic solution. Additionally, when the bubble is confined by the high pressure from the ambient medium, we classify it as the cooling, pressure-confined (CP) mode.
• 6.  
  No wind (NW) mode: The mass-loss rate is extremely low (e.g., 10⁻⁴ M_⊙ yr⁻¹) and the ambient medium is very dense (e.g., n = 10³ cm⁻³ here), so freely expanding supersonic wind is completely inhibited. The wind solution does not meet the outflow pressure criterion of Cantó et al. (2000) (see below).
• 7.  
  Momentum-conserving (MC) mode: The mass-loss rate is high and the velocity is low, causing catastrophic cooling so that the wind is suppressed. We classify as MC mode those models having a mean initial temperature in region (a) that is a factor of 3 below the expected mean adiabatic temperature within a distance of 1 pc from the cluster boundary.

The classifications fall into three categories: adiabatic and quasi-adiabatic (AW, AB, AP), radiatively cooling (CC, CB, CP), and suppressed superwinds (NW, MC). The suppressed superwinds are the extreme cases. The NW mode represents the case where the mass-loss rate is too low relative to the ambient density and a wind cannot launch. As mentioned by Cantó et al. (2000), the existence of the supersonic wind depends on the following outflow pressure:

$$\begin{eqnarray}&&{P}_{\mathrm{out}}\equiv \displaystyle \frac{1}{\gamma }{\left(\displaystyle \frac{\gamma -1}{\gamma +1}\right)}^{\tfrac{1}{2}}{\left(\displaystyle \frac{\gamma +1}{2}\right)}^{\tfrac{\gamma }{\gamma -1}}\displaystyle \frac{{\dot{M}}_{\mathrm{sc}}{V}_{\infty }}{4\pi {R}_{\mathrm{sc}}^{2}},\end{eqnarray} \tag{ 14 }$$

To have the supersonic solution, it is necessary to have the outflow pressure P_out higher than the ambient thermal pressure P_amb in the surrounding ambient medium. In the ideal gas condition, the ambient thermal pressure P_amb is a function of the ambient temperature and density as P_amb = k_B T_amb ρ_amb/μ, where k_B is Boltzmann's constant and μ is the mean mass per particle of the gas. This provides a physical explanation for the dependence of the wind existence on the mass-loss rate, the wind terminal velocity, and the ambient density, for a given cluster radius (R_sc = 1 pc) and ambient temperature. Outflow with very low mass-loss rate and/or terminal velocity does not support the supersonic solution, so there is no freely expanding supersonic wind. From Equation (14), it can be seen that supersonic conditions also depend on the cluster radius (R_sc), so at a very large value of R_sc, depending on ${\dot{M}}_{\mathrm{sc}}$ and V_∞, the existence of freely expanding supersonic winds is also impossible.

Superwinds are also suppressed by the inverse situation where the mass-loss rate is large and the velocity low, causing strong catastrophic cooling. The MC mode corresponds to conditions like these, where there is no significant thermal energy contributing to the outflow, and it roughly proceeds in a momentum-conserving mode only (see, e.g., Ostriker & McKee 1988; Koo & McKee 1992a, 1992b). For such models with high mass-loss rates and low wind speeds (e.g., V_∞ = 250 km s⁻¹ for 10⁻² M_⊙ yr⁻¹), it might often be the case that a supersonic wind cannot be launched at the injection radius R_sc (e.g., Silich & Tenorio-Tagle 2018). This could arise in situations such as those studied by Silich & Tenorio-Tagle (2018), where hot shocked winds from individual stars strongly cool down inside the cluster before a superwind can be launched into the ambient medium. Previously, it was shown that high mass loading can contribute to strong radiative cooling in the hot shocked wind (Lochhaas et al. 2018; Silich et al. 2020), which may also suppress or delay the development of a star cluster superwind (Silich et al. 2020).

The parameter space for the different regimes can be seen in Figure 4. This shows the temperature discrepancy factor (f_T ≡ T_w/T_adi) defined as the mean value of the expanding wind temperature T_w over the predicted analytic adiabatic temperature T_adi ∝ r^−4/3 for the freely expanding supersonic wind, region (a). A value of f_T lower than 0.75 is associated with strongly radiative cooling (CC and CB modes). A value for the temperature discrepancy factor of ∼1 generally corresponds to the adiabatic cooling mode with 0.75 < f_T < 1.25 (Model AB). It can be seen that, as expected, an increase in the mass-loss rate decreases the temperature discrepancy factor, because the higher density increases cooling. Increasing the metallicity also reduces f_T , thus also increasing radiative cooling (e.g., Silich et al. 2004). The MC models represent the extreme case, where there is no thermal energy contributing to the outflow. We also see that for a given ${\dot{M}}_{\mathrm{sc}}$, catastrophic cooling occurs for lower outflow velocities. While strong cooling is also promoted by higher metallicity, Figure 4 shows that ${\dot{M}}_{\mathrm{sc}}$ and V_∞ are more important over the modeled parameter space.

Zoom In Zoom Out Reset image size

The pressure-confined models (AP) tend to have higher temperatures, as expected since the adiabatic expansion is reduced, therefore elevating the temperatures in the interior regions according to the ideal gas law. Two models have f_T > 1.25. They are identified in Figure 4 as ${\dot{M}}_{\mathrm{sc}}={10}^{-4}$ M_⊙ yr⁻¹, V_∞ = 1000 km s⁻¹, n_amb = 1000 cm⁻³, Z/Z_⊙ = 1 and ${\dot{M}}_{\mathrm{sc}}=0.224\times {10}^{-3}$ M_⊙ yr⁻¹, V_∞ = 736 km s⁻¹, n_amb = 1000 cm⁻³, Z/Z_⊙ = 0.125. These are models at the high-density limit, where bubble growth is inhibited the most.

As shown in Figure 5, the density discrepancy factor (f_n ≡ n_w/n_adi), defined as the mean value of the expanding wind density n_w over the adiabatic density n_adi, is computed for the freely expanding supersonic wind (region (a)). It can be seen that the density profiles mostly follow the adiabatic solution. However, the regions with strong cooling (CC, CB) have slightly enhanced densities, as expected. High densities could also lead to high ambient pressures, which also lead to slower bubble growth, sometimes resulting in the adiabatic, pressure-confined (AP) models.

Zoom In Zoom Out Reset image size

We also see that the wind mechanical power $\tfrac{1}{2}{\dot{M}}_{\mathrm{sc}}{V}_{\infty }^{2}$ does not itself determine the outflow's status with respect to catastrophic cooling; some models with the weakest wind power (upper left group in Figure 4) are still adiabatic. Rather, ${\dot{M}}_{\mathrm{sc}}$ and V_∞ are independently the critical parameters in establishing strong cooling. Moreover, we see that there are more catastrophic cooling models that are still capable of supporting hot bubbles (CB) than those without (CC). Thus, the presence of a hot, X-ray superbubble does not necessarily indicate adiabatic feedback.

Figure 4 shows that for the models with fiducial metallicity scalings for the effective mass-loss rate and wind velocity (${\dot{M}}_{\mathrm{sc}}={10}^{-2}$ M_⊙ yr⁻¹ and V_∞ = 1000 km s⁻¹), essentially all the models show conventional adiabatic feedback. Reducing ${\dot{M}}_{\mathrm{sc}}$ does not change this scenario until the outflow simply cannot penetrate the ambient density. However, reducing the wind velocity, or otherwise reducing the kinetic heating efficiency, destabilizes these models to catastrophic cooling. Reducing the wind density can mitigate this effect, but only to a point, as seen in these model grids. In practice, regions that tend to be radiation dominated are probably more likely to be mass-loaded owing to photoevaporation. A smaller outflow region, which keeps the gas close to the SSC, also promotes these effects. On the other hand, if catastrophic cooling takes place inside the launch radius (e.g., Silich et al. 2020), it could result in a very small effective ${\dot{M}}_{\mathrm{sc}}\sim {10}^{-4}\ {M}_{\odot }{{yr}}^{-1}$. This could potentially generate the parameter space of models in that regime, which are otherwise unrealistic for stellar wind outflows.

In this section, we present emission-line emissivities and luminosities determined using CLOUDY, for the two different models in Section 4, namely, PI (pure photoionization) and CPI (photoionization + hydrodynamic collisional ionization).

Figure 6 (top panels) shows the emissivities calculated for various emission lines as a function of radius for our photoionization models with the two different ionization schemes (PI and CPI; right panels). The PI emissivity profiles, which were made using the Starburst99 SED and luminosity applied to the neutral density profile from the hydrodynamic simulation, are shown in the top panel. The bottom panel presents the emissivities from the CPI model built using both the Starburst99 SED ionizing source and hydrodynamic ionization. Similarly, the ionic fractions for our photoionization models with these two different (PI and CPI) models are presented in Figure 6 (bottom panels).

The CPI models are used to calculate line emissivities in what follows. We classify our photoionization models according to the nebular (H⁺) optical depth of the swept-up shell. Optically thick models are shown with boldface in Figure 4. Models are optically thin if the ambient medium beyond the shell (region (c) in Figure 2) is photoionized to produce H⁺ by the SSC. We see that the optically thick models are those in the densest (n_amb ≳ 10²) ambient media, for which dense, optically thick swept-up shells develop more quickly. However, we caution that our 1D models do not account for shell clumping, which can greatly alter the escape fraction of ionizing radiation. We also caution that MAIHEM does not implement radiative transfer, which may be a significant effect in high-density models. Thus, the exact boundary of the ionized edge is approximate, and the ionization structure of the interface region between the optically thin and thick conditions is also approximated. Thus, for the optically thick CLOUDY models, we adopt the radial temperature profile from the PI model from the radius at which the shell becomes isothermal according to the MAIHEM simulation, roughly 1–2 pc outward from the inner boundary of the dense shell (see Section 3).

Observations of distant, luminous complexes may be biased toward the central, high surface brightness regions, thereby approximating density-bounded observations transverse to the line of sight. To allow for the calculation of partially and fully radiation-bounded models, the total luminosity L_λ of each emission line at wavelength λ is calculated from its volume emissivity _λ (r) as follows (see Appendix A):

$$\begin{eqnarray}&&{L}_{\lambda }=4\pi {\int }_{R=0}^{{R}_{\mathrm{aper}}}\left[{\int }_{r=R}^{{R}_{\max }}\displaystyle \frac{{\epsilon }_{\lambda }(r)}{\sqrt{{r}^{2}-{R}^{2}}}{rdr}\right]{RdR},\end{eqnarray} \tag{ 15 }$$

where r is the radial distance from the center, R the projected distance from the center, R_aper the radius of the circular projected aperture within the object observed, and ${R}_{\max }$ the maximum radius of the object in the line of sight. Thus, for a fully "density-bounded" model, we set ${R}_{\max }={R}_{\mathrm{shell}}$, i.e., the maximum radius of the dense shell (region (c) in Figure 2), and R_aper = R_shell. For a fully "radiation-bounded" H ii region, ${R}_{\max }$ and R_aper are set to the Strömgren radius ${R}_{\mathrm{Str}}$ predicted by CLOUDY. Setting R_aper = R_shell and ${R}_{\max }={R}_{\mathrm{Str}}$ is radiation bounded in the line of sight, but otherwise density bounded. This geometry is described as a "partially density-bounded" model in what follows. Appendix A fully describes the calculations based on this geometry.

Table 3 partially lists the total luminosities derived from emissivities calculated by CLOUDY for radiation-bounded, partially density-bounded, and density-bounded models with the two different cases (PI and CPI). The tables for all the photoionization calculations are presented in the machine-readable format in Appendix B.

5.1. Diagnostic Diagrams

Here we explore the emission-line parameter space generated by our CPI models in a few diagnostic diagrams, including the locus for catastrophic cooling models. We use our total luminosities derived from CLOUDY emissivities to produce diagnostic diagrams for different metallicities, wind terminal velocities, ambient densities, and mass-loss rates.

Figure 7 shows "BPT diagrams" (Baldwin et al. 1981) based on [O iii] λ5007/Hβ versus [S ii] λ λ6717, 6731/Hα and [N ii] λ6584/Hα for our fully radiation-bounded models with optically thin and thick shells plotted by filled and open circles, respectively. As noted above, given the uncertainties when the transition between optically thin and thick conditions takes place in the shell, the line emissivities for optically thick models are somewhat uncertain. These diagrams are widely used to distinguish between active galactic nuclei (AGNs) and starburst galaxies (e.g., Veilleux & Osterbrock 1987; Osterbrock et al. 1992; Kewley et al. 2001, 2006, 2013; Groves et al.2004, 2006). Similarly, Figure 8 shows the same models as in Figure 7 but with the emission-line luminosities calculated as density bounded at the shell outer radius transverse to the line of sight, as described above.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In general, we see that both sets of models in Figure 7 and 8 show line ratios typical of photoionized H ii regions. This applies to both adiabatic and strongly cooling conditions. The models are highly excited, although slightly below the starburst boundaries defined by Kewley et al. (2001) and Kauffmann et al. (2003). This is likely caused by the fact that the ionization parameter is diluted by the outflows in our models, which displace the gas to larger radii at the shells. Higher-density models also have higher ionization parameters, as evidenced by their low [N ii] and [S ii] emission. As expected, at higher metallicity, [O iii]/Hβ decreases owing to stronger cooling. The partially density-bounded models in Figure 8 are similar to the fully radiation-bounded ones in Figure 7, but lacking the strongest [N ii] and [S ii] emission, which ordinarily dominates in the outer regions. Since the wind itself is more highly ionized and low density, the contributions from kinematic heating are not apparent in the species shown by the BPT diagrams.

Similar to the optical BPT diagrams, UV diagnostics diagrams based on ultraviolet emission lines such as C iv λ λ1549, 1551, He ii λ1640, O iii] λ λ1661, 1666, C iii] λ λ1907, 1909, He ii λ1640 have been employed to distinguish between AGNs and starburst galaxies (Allen et al. 1998; Feltre et al. 2016; Gutkin et al. 2016; Hirschmann et al. 2019) and radiative shock regions (Allen et al. 1998; McDonald et al. 2015; Hirschmann et al. 2019). We select sets of those UV diagnostic line ratios that are typically reported in star-forming regions to make UV diagnostic diagrams. Fully radiation-bounded models are presented in Figures 9 and 11, and models that are density-bounded transverse to the line of sight are shown in Figures 10 and 12.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 presents the UV diagnostic diagrams for O iii] λ λ1661, 1666/He ii λ1640 versus C iv λ λ1548, 1551/C iii] λ λ1907, 1909 and C iv λ λ1548, 1551/He ii λ1640, for our fully radiation-bounded models. In general, these line ratios increase for decreasing metallicity, due to the increasing nebular electron temperature in metal-poor systems. However, at Z/Z_⊙ = 1, strong adiabatic feedback generates hot bubbles with strong C iv emission in the interface with the dense shell, displacing some models for V_∞ = 1000 km s⁻¹ to higher C iv ratios relative to other models. For all other models, the wind heating is unimportant for these emission-line ratios. This highlights the weakness of stellar winds at low metallicity.

These effects are more pronounced in Figure 10, which is the same as Figure 9, but for models that are partially density bounded. Since these models are weighted toward emission in and near the shell, the C iv emissivity is enhanced in models with significant adiabatic feedback, such as those at lower density and higher velocity. On the other hand, O iii] is more sensitive to the overall electron temperature, which decreases with increasing ambient density.

The highly ionized doublet O vi λ λ1032, 1038 originates from shocks and hydrodynamic collisional ionization, rather than photoionization. Gray et al. (2019b) suggested that it could be a useful diagnostic of catastrophic cooling flows. In Figures 11 and 12, we examine O vi λ λ1032, 1038/He ii λ1640 versus C iv λ λ1549, 1551/He ii λ1640, for fully radiation-bounded and partially density-bounded models transverse to the line of sight, respectively. In our models, the O vi emission is produced primarily in the central region of the free-expanding wind, nearest to the SSC, and in the dense shell interfaces to the hot bubble region (see Figure 6, top panels).

Zoom In Zoom Out Reset image size

Figures 11 and 12 show that denser winds increase the O vi emission through strong cooling, which correlates strongly with ${\dot{M}}_{\mathrm{sc}}$, as well as showing a slight anticorrelation with V_∞. These effects are enhanced in the partially density-bounded plots in Figure 12. This is consistent with the suggestion by Paper I that O vi is associated with catastrophic cooling, but we also see that adiabatic models also show significant O vi. Thus, while O vi is indeed enhanced in most catastrophic cooling models, it is difficult to use as a definitive diagnostic of these conditions without other constraints on the parameter space. Indeed, for the most strongly cooling models, in which the hot bubble region is eliminated, the temperature is too low to generate O vi at all.

Zoom In Zoom Out Reset image size

We caution that, as demonstrated by Gray et al. (2019b), the O vi emission line can be strong in the outflow owing to NEI states. Moreover, Paper I shows that NEI conditions are expected to contribute to more highly ionized states than in CIE, especially in strongly cooling flows. This is also seen, for example, in SN remnants (e.g., Patnaude et al. 2009, 2010; Zhang et al. 2019).

In Figures 7–12, we also see that optically thick models typically appear in denser ambient media. Figure 7 shows that the ambient density of n_amb = 10³ cm⁻³ results in weaker singly-ionized emission lines [N ii] and [S ii]. As seen in Figure 9, an increase in the metallicity (Z/Z_⊙) mostly tends to decrease the UV line ratios O iii]/He ii and C iv/He ii, since the nebular temperature decreases. It is also not clear that enhanced C iv emission is strongly associated with catastrophic cooling models as suggested in Paper I. However, future work with NEI models is needed to clarify this.

As noted above, our CLOUDY models include dust grains with M_d/M_Z = 0.2 found by De Vis et al. (2019) for evolved galaxies. Dust grains could also form in shocked ejecta from SN explosions (e.g., Nozawa et al. 2010) in starburst-driven superwinds. On the other hand, the dust-to-metal ratio could be lower for unevolved galaxies. Although the optical diagnostic line ratios plotted in the BPT diagrams are not greatly affected by the presence of dust grains, they can have significant changes on the UV line ratios, since dust grains slightly decrease radiation fields. Moreover, metal depletion onto dust grains reduces the ratios of C iv λ λ1549, 1551/He ii λ1640 and O vi λ λ1032, 1038/He ii λ1640. Thus, dust grains and depletion onto them are additional sources of uncertainties in these line ratios, as well as in NEI cooling functions (e.g., Richings et al. 2014).

We compare our results with optical and UV observations of nearby and distant starburst galaxies, which were found to have detected observations of He ii λ1640, C iv λ1550, O iii] λ1664, and C iii] λ1908. The observations of low-metallicity compact dwarf galaxies studied by Berg et al. (2016) have mean oxygen abundance of $12+\mathrm{log}($ O/H ) ≈ 7.6 ($\hat{Z}\approx 0.125$) and electron temperatures T_e([N ii]) ≈ 15,000 K. Nearby starburst dwarf galaxies studied by Senchyna et al. (2017) have similar parameters, but with $12+\mathrm{log}($O/H) ≲ 8.3 ($\hat{Z}\lesssim 0.6$). More distant starburst galaxies have been extensively analyzed using optical nebular lines, and some have been recently studied through rest-frame UV spectroscopic observations (Stark et al. 2014; Patrício et al. 2016; Steidel et al. 2016; Amorín et al. 2017; Vanzella et al. 2016; Saxena et al. 2020). In particular, Amorín et al. (2017) considered 10 starburst galaxies at redshift 2.4 < z < 3.5 with mean metallicity of $12+\mathrm{log}($ O/H ) ≈ 7.6, for which they measured Lyα λ1216, C iv λ1550, He ii λ1640, O iii] λ1664, and C iii] λ1908. These objects have similar properties and line ratios to the samples of Berg et al. (2016) and Senchyna et al. (2017), and they are also included in Figures 9 and 10. Their properties are largely the same as those of the local sample for the considered line ratios. The galaxy samples occupy the parameter space used for our models, as shown in Figures 7 to 10.

It can be seen in Figures 7 and 8 that, for the appropriate metallicity within Z/Z_⊙ = 0.125 and 0.5, the data agree more with our fully radiation-bounded models than with the partially density-bounded models. This implies that the conditions in these objects do not require significant density bounding for the integrated nebular emission. However, since these apertures are on the order of a few kiloparsecs, we note that the observed singly-ionized emission lines could be due to the diffuse ISM, and so some density bounding is not ruled out. These samples have similarities to the Green Pea galaxies, which tend to have higher ionization parameters, with weaker emission from the low-ionization species in Figures 7 and 8.

Figures 9–10 also show the UV observational data from Berg et al. (2016) and Senchyna et al. (2017) plotted on the UV diagnostic diagrams. We again see reasonably good correspondence between the data and models for Z/Z_⊙ = 0.5 and 0.25. In Figure 9, the models appear to slightly overpredict the O iii]/He ii ratio. For these radiation-bounded conditions, the discrepancy may be due to the CIE approximation in our models. As shown in Paper I, NEI conditions for such outflows generate slightly higher ionization states, which would reduce O iii]/He ii and C iii]/C iv. Furthermore, Figure 10 also shows that only a small degree of density bounding for the high-density models can alleviate the discrepancy. This density bounding could be an observational bias, as suggested earlier, or it could be at least partially real, since similar galaxies have been found to be Lyman continuum (LyC) emitters (e.g., Izotov et al. 2016, 2018).

The O vi λ λ1032, 1038 emission and absorption lines are more readily associated with active galaxies and AGN feedback, rather than starbursts (e.g., Baldwin et al. 2003; Hainline et al. 2011; Tripp et al. 2008; Savage et al. 2014; Bielby et al. 2019). We do not typically expect O vi emission in ordinary H ii regions, since these highly ionized lines originate from a hot gas phase with T ∼ 3 × 10⁵ K. Generally in starbursts, this species is only seen in absorption, associated with adiabatic outflows (e.g., Heckman et al. 2001). But recently, Hayes et al. (2016) spatially resolved O vi λ λ1032, 1038 emission in a nearby, intense starburst galaxy (SDSS J115630.63+500822.1, hereafter J1156) through Hubble Space Telescope (HST) imaging, having spectroscopically observed this emission with the Far Ultraviolet Spectroscopic Explorer (FUSE). Moreover, they also detected a blueshifted O vi absorption line associated with an average outflow velocity of 380 km s⁻¹. From the O vi column density, they found an O vi mass of 3 × 10⁴ M_⊙. The starburst properties of this galaxy are similar to those of our comparison data sets, and it has metallicity $7.87\lt 12+\mathrm{log}($O/H) < 8.34 (i.e., $0.23\lt \hat{Z}\lt 0.7$).

Chisholm et al. (2018) confirmed that photoionization models of galactic outflows with CLOUDY cannot reproduce the observed O vi absorption in a high-redshift starburst, so they again point to the standard model that O vi originates from a mass-loading interface between the photoionized region and unobserved hot wind. However, Gray et al. (2019b) found that NEI models can produce a range of O vi column densities in dense outflows themselves, confirming that relatively high O vi emission can also be due to hydrodynamic heating. Moreover, Paper I also reproduced strong O vi and C iv emission from strongly cooling winds, so they suggested that these lines could be used to trace catastrophically cooling superwinds in starburst regions. In Figure 11, we see that the O vi emission is indeed enhanced by increasing the mass-loss rate and metallicity, and it is typically stronger in catastrophic cooling models. However, the O vi emission plummets in the models at Z/Z_⊙ = 1 and ${\dot{M}}_{\mathrm{sc}}$=10⁻³ M_⊙ yr^–1, where strong radiative cooling occurs. Moreover, O vi is also seen in adiabatic models, so additional constraints on the parameter space are necessary when using this line to evaluate the presence of catastrophic cooling.

The local extreme emission-line galaxies are regarded as analogs of higher-redshift starbursts. Various samples of such objects have been reported, and they have similar properties and emission-line ratios to the comparison samples above from Berg et al. (2016) and Senchyna et al. (2017). For example, the recent deep VANDELS spectroscopic survey of He ii emitting galaxies by Saxena et al. (2020) identified eight objects with high enough signal-to-noise ratio at z ∼ 2.5–3.5 to detect He ii λ1640, O iii] λ λ1661, 1666 lines, and C iii] λ1909 in emission. The metallicities of these He ii emitters are estimated to be in the range Z/Z_⊙ = 0.125–0.25, similar to the low-redshift samples. Other UV observations of intermediate-redshift sources have similar physical conditions. Deep spectroscopic observations of low-mass, gravitationally lensed galaxies at z ∼ 1.5–3 conducted by Stark et al. (2014) identified moderately metal-poor gaseous environments with ∼0.2–0.5 Z_⊙, and MUSE integral field spectroscopy of a young lensed galaxy at z = 3.5 (SMACS J2031.8–4036) by Patrício et al. (2016) revealed similar physical conditions (T_e ∼ 15,600 K) and metallicity ( ≲ 0.16 Z_⊙).

These highly excited starbursts appear to be linked to LyC emission. The extreme Green Pea galaxies are to date the only class of objects known to be confirmed LyC emitters (e.g., Izotov et al. 2016, 2018), and they are characterized by their extreme nebular excitation (Cardamone et al. 2009). At intermediate redshift, de Barros et al. (2016) investigated the UV spectrum of an LyC emitter candidate at z = 3.2, called Ion2, which has $12+\mathrm{log}($O/H) ≈ 8.1 ($\hat{Z}\approx 0.4$) and strong Lyα [O iii] λ λ4959, 5007 and C iii] λ λ1907, 1909. Ion2 exhibits an extremely high O iii]/[O ii] ratio > 10, which is much higher than expected according to the SFR–M_⋆–Z dependence of O iii]/[O ii] ratios derived by Nakajima & Ouchi (2014). De Barros et al. (2016) suggest that this could be due to a density-bounded scenario. Indeed, primeval galaxies at the time of cosmic reionization (z > 7) also have these extreme emission-line properties, and similarly are very compact, with high star formation rates and lower masses (e.g., Grazian et al. 2015; Shibuya et al. 2015; Tasca et al. 2015). Stark et al. (2015) analyzed the UV observations of a gravitationally lensed, Lyα emitter at z = 7.045 (A1703-zd6) and having $12+\mathrm{log}($O/H) ≈ 7.04. They detected C iv λ1548 emission with FWHM ≲ 125 km s⁻¹. Although this is similar to measurements in narrow-lined AGNs at lower redshift, it could also originate from a starburst region powered by young, very hot, metal-poor stars.

Catastrophic cooling conditions have been suggested to be present in some of these extreme starbursts (e.g., Silich et al. 2020; Jaskot et al. 2019; Oey et al. 2017) based on kinematic and other evidence. If this is indeed the case, then our results would imply that most likely the heating efficiency in these systems is significantly reduced relative to the fiducial scalings for stellar wind velocities (Section 4). It is also possible that mass loading plays a role, but the degree of mass loading required to induce catastrophic cooling is high. Since these objects tend to be metal-poor, it is unlikely that metallicity plays a significant role in promoting cooling. If this is indeed the case, then our results in Section 4 would imply that either (1) the heating efficiency in these systems is low, (2) mass loading is high, or (3) the winds are suppressed within the injection radius R_sc. The first option would imply that the effective wind velocity is reduced by at least a factor of 2–4, since stellar winds are on the order of 1000–2000 km s⁻¹, even at low metallicity. This would correspond to a decrease in kinetic heating efficiency of 0.25–0.06 and could be caused by a variety of factors like entrainment of clouds, turbulence, and magnetic fields. The second option would similarly require an increase in effective ${\dot{M}}_{\mathrm{sc}}$ by an order of magnitude, which could be accomplished by mass loading. The third option would require both weak velocities and high ${\dot{M}}_{\mathrm{sc}}$ in order to suppress the wind within R_sc. This has been suggested to be the case in NGC 5253-D1 by Silich et al. (2020).

Thus, our modeled emission-line spectra and outflow parameter space may be relevant to LyC emitters and primordial galaxies. Our comparison with the local extreme starbursts shows that our models for Z/Z_⊙ = 0.5 and 0.25 are reasonably good at matching the observed emission-line data for C iv, C iii], O iii], and He ii when invoking a slight density bounding at high ambient densities, on the order of 100 cm⁻³. This density bounding would be expected for LyC emitters, although it could also simply be an observational bias favoring the brightest, swept-up gas. Or, as discussed in Section 6, the objects are fully radiation bounded but with NEI conditions in the outflow that cause the modest variation from our modeled emission-line ratios. Both scenarios imply a significant contribution from kinematic heating to generate the observed high ionization states. Indeed, catastrophic cooling has been suggested to be linked to LyC-leaking galaxies (Jaskot et al. 2017, 2019).

As suggested by Paper I, catastrophic cooling models as a group tend to show higher emission in C iv and O vi, especially for the CB models. But for stronger cooling in the CC models, the temperature is too low to support strong emission in these ions. We also see that adiabatic models also produce C iv and O vi, complicating their use as diagnostics of catastrophic cooling, although additional work is needed to evaluate this. Since our models are still limited in parameter space, our comparisons with observations are not conclusive but illustrate the feasibility of the catastrophic cooling and density bounding as important factors in observed line emissivities. Further study of known LyC emitters with resolved nebular data, such as Haro 11 (e.g., Micheva et al. 2020; Keenan et al. 2017), is also needed to develop robust nebular diagnostics of these conditions.

In this paper, we employ the nonequilibrium atomic chemistry and cooling package MAIHEM (Gray et al. 2015; Gray & Scannapieco 2016; Gray et al. 2019b) to investigate how superwinds driven by SSCs are subject to radiative cooling. We assume an SSC with mass 2.05 × 10⁶ M_⊙, radius of 1 pc, and age of 1 Myr, and we model a range of ambient density (n_amb = 1,...,10³ cm⁻³). At solar metallicity, the superwind is parameterized by mass loss ${\dot{M}}_{\mathrm{sc}}$ = 10⁻¹,...,10⁻⁴ M_⊙ yr^–1 and wind terminal velocity V_∞ = 250, 500, 1000 km s⁻¹. At lower metallicity (Z/Z_⊙ = 0.5, 0.25, 0.125), ${\dot{M}}_{\mathrm{sc}}$ and V_∞ are scaled by the metallicity dependencies, i.e., ${\dot{M}}_{\mathrm{sc}}\propto {Z}^{0.72}$ and V_∞ ∝ Z^0.13. The superwind model is implemented following Gray et al. (2019a), according to the radiative solution presented by Silich et al. (2004).

Our resulting grid of superwind models demonstrates where strongly radiative cooling occurs in this parameter space. The superwind structures produced by our hydrodynamic simulations are classified according to their departure from the expected adiabatic temperature, as well as the formation of the characteristic bubble within the simulated time frame. Our models are parameterized in terms of ${\dot{M}}_{\mathrm{sc}}$ and V_∞, which account for mass loading and/or heating efficiency effects.

We find that decreasing V_∞ or, equivalently, the kinetic heating efficiency significantly enhances radiative cooling effects for a given fixed ${\dot{M}}_{\mathrm{sc}}$ and n_amb (Figure 4). Similarly, high mass-loss rates (${\dot{M}}_{\mathrm{sc}}$) or, equivalently, high mass loading can lead to catastrophic cooling in slower stellar winds (lower V_∞). Where both of these conditions apply, radiative cooling strongly dominates and prevents any dynamically significant thermal energy contribution, resulting in either momentum-conserving outflows or winds failing to launch from the cluster injection radius R_sc. Increasing Z also promotes cooling, but it is not as strong an effect as the wind parameters.

We also find that the presence of the distinctive, hot bubble morphology is not always a reliable indicator of adiabatic status for the outflow. The hot bubble accompanies both adiabatic outflows and ones with strong cooling flows. Therefore, the existence of a hot superbubble is not necessarily indicative of adiabatic feedback. Similarly, some adiabatic models take time to develop the bubble morphology. Thus, for young systems with ages ≲1 Myr, the lack of a hot, X-ray bubble is not necessarily an indicator that the system is not adiabatic.

We calculate line emissivities for the density and temperature profiles produced by our superwind models using the photoionization code CLOUDY (Ferland et al. 1998, 2013, 2017). The line luminosities are computed from the volume emissivities for fully radiation-bounded models, and models that are density-bounded transverse to the line of sight. We construct a few optical and UV diagnostic diagrams from the predicted line luminosities (Figures 7–12) to compare with observed emission lines from extreme starbursts. As noted above, dust is also a factor in evaluating the UV line ratios.

As suggested by Paper I, O vi and C iv are stronger in modest catastrophic cooling flows, especially the former. For extreme cooling, the temperature is too low to generate strong emission in these ions. Since adiabatic models also produce these emission lines, their use as diagnostics of catastrophic cooling is complicated and requires careful analysis of the parameter space. Further study is needed to identify unambiguous diagnostics.

Our optical and UV diagnostic line ratios predicted by radiation-bounded models are in general agreement with observations of nearby and distant star-forming galaxies having similar physical properties (see Figures 7–12). The modest discrepancies observed in our UV diagnostic diagrams in C iv, C iii], O iii], and He ii are consistent with minor density bounding, implying LyC escape, as has been suggested for similar galaxy populations. Or, the observed density bounding could be an observational bias toward the dense, piled-up shells. NEI conditions may also be responsible. However, all of these interpretations imply significant emission from kinetically heated thermal components in our modeled outflows.

We caution that our 1D hydrodynamic simulations do not reproduce thermal and dynamical instabilities that create clumping in the free-expanding superwind and shell. This would significantly affect the outflow morphology and evolution, and, by extension, the expected emission-line spectra. Instability-induced gas clumping within the outflow may be important (e.g., Jaskot et al. 2019) and needs to be investigated through 2D hydrodynamic simulations. Other important parameters remain to be explored, including the effect of the SSC size, ambient density gradient, system age, and NEI. Our photoionization calculations are implemented using the typical dust-to-metal mass ratio associated with evolved galaxies, whose metallicity dependence is not well understood. Future computations are required to investigate these effects in order to better understand the occurrence and properties of catastrophic cooling in starburst regions.

We are grateful to Sergiy Silich for useful discussions and comments on the manuscript. We thank the referee for comprehensive comments and suggestions. M.S.O. acknowledges NASA grants HST-GO-14080.002-A and HST-GO-15088.001-A.

The hydrodynamics code FLASH used in this work was developed in part by the DOE NNSA ASC- and DOE Office of Science ASCR-supported Flash Center for Computational Science at the University of Chicago. Analysis and visualization of the FLASH simulation data were performed using the yt analysis package (Turk et al. 2011).

Software: FLASH (Fryxell et al. 2000), yt (Turk et al. 2011), CLOUDY (Ferland et al. 2017), Starburst99 (Leitherer et al. 2014), NumPy (Harris et al. 2020), SciPy (Virtanen et al. 2020), Matplotlib (Hunter 2007), HDF5 (Folk et al. 2011), hypre (Falgout & Yang 2002).

For comparison of CLOUDY results with observations, we need to generate the projected 2D surface brightness and the total luminosity from the volume emissivity created by CLOUDY for each emission line. We can then compare the observed spatially resolved flux with the projected 2D surface brightness and the observed integrated flux with the total luminosity.

We set the program CLOUDY to generate the volume emissivity _λ (r) at a given radius r for each emission line at wavelength λ. Figure 13 shows a schematic of the transformation from the volume emissivity _λ (r) to the surface brightness I_λ (R) projected onto the x-y plane perpendicular to the line of sight along the z-axis. We assume that the H ii region is optically thin and spherically symmetric. To calculate the surface brightness of a spherically symmetric object at each point on the x and y projected plane, we employ the SciPy package to integrate the volume emissivity _λ (r) over the z-axis, ${I}_{\lambda }(R)=2{\int }_{z\,=\,0}^{{z}_{\max }}{\epsilon }_{\lambda }(r){dz}$, where $z=\sqrt{{r}^{2}-{R}^{2}}$ and ${dz}={rdr}/\sqrt{{r}^{2}-{R}^{2}}$, so the projected surface brightness at the projected radius from the center of the object is calculated by (see, e.g., Sarazin 1986; Cimatti et al. 2019)

$$\begin{eqnarray}&&{I}_{\lambda }(R)=2{\int }_{r=R}^{{R}_{\max }}\displaystyle \frac{{\epsilon }_{\lambda }(r)}{\sqrt{{r}^{2}-{R}^{2}}}{rdr},\end{eqnarray} \tag{ A1 }$$

where r is the radial distance of the volume emissivity from the center of the object, R is the projected radius of a point on the 2D projected plane from the center of the object, and ${R}_{\max }$ is the maximum radius of the geometry where the integral is performed over and the volume emissivity _λ (r) should also be available at the maximum radius ($r={R}_{\max }$).

Zoom In Zoom Out Reset image size

We note that Equation (A1) is an Abel integral, so it can be inverted and deprojected to the volume emissivity as a function of radius as follows (see, e.g., Cavaliere 1980; Longair 2008):

$$\begin{eqnarray}&&{\epsilon }_{\lambda }(r)=-\displaystyle \frac{1}{\pi r}\displaystyle \frac{d}{{dr}}{\int }_{R=r}^{{R}_{\max }}\displaystyle \frac{{I}_{\lambda }(R)}{\sqrt{{R}^{2}-{r}^{2}}}{RdR}.\end{eqnarray} \tag{ A2 }$$

The total luminosity integrated over a circular aperture entirely or partially covering a spherical symmetric H ii region (see Figure 13) with the surface brightness I_λ (R) projected on the 2D plane is

$$\begin{eqnarray}&&{L}_{\lambda }={\int }_{\varphi =0}^{2\pi }{\int }_{R=0}^{{R}_{\mathrm{aper}}}{I}_{\lambda }(R){RdRd}\varphi ,\end{eqnarray} \tag{ A3 }$$

where R_aper is the radius of the circular aperture where the integral is performed over. Substituting Equation (A1) into the above integral (A3), we restore Equation (15).

We use Equation (A1) to construct the 2D projected surface brightness that is comparable with the observed spatially resolved emission-line flux. Moreover, we utilize Equation (A3) to calculate the total luminosity, which can be compared with the observed integrated emission-line flux.

To model the observed emission that is radiation bounded in the line of sight but density bounded in transverse directions, we set R_aper = R_shell and ${R}_{\max }={R}_{\mathrm{Str}}$, i.e., a partially density-bounded model; however, ${R}_{\mathrm{aper}}={R}_{\max }={R}_{\mathrm{shell}}$ for a fully density-bounded model, and ${R}_{\mathrm{aper}}={R}_{\max }={R}_{\mathrm{Str}}$ for fully radiation-bounded model, where R_shell is the outer radius of the dense shell and ${R}_{\mathrm{Str}}$ is the Strömgren radius (see Figure 13).

The compressed (tar.gz) files of the interactive figure (64 images) of Figure 2, the interactive animation (64 videos) of Figure 3, and the interactive figure (192 images) of Figure 6 are available in the electronic edition of this article and are archived on Zenodo (doi:10.5281/zenodo.4989577). These interactive figures and animation are also hosted on this URL. ⁴

The compressed (tar.gz) file containing the six machine-readable tables for the emission-line data (partially presented in Table 3) is available in the electronic edition of this article. Each file is named as table_case_bound.dat such as table_CPI_radi.dat, where case is for the ionization case (PI: purely photoionization; CPI: photoionization and hydrodynamic collisional ionization) and bound for the optical depth model (radi: fully radiation bounded; pden: partially density bounded; dens: fully density bounded). Each file contains the following information:

• 1.  
  metal: metallicity $\hat{Z}\equiv Z/$ Z_⊙ = 1, 0.5, 0.25, and 0.125.
• 2.  
  dMdt: mass-loss rate ${\dot{M}}_{\mathrm{sc}}={10}^{-1}$, 10⁻², 10⁻³, and ${10}^{-4}\times {\hat{Z}}^{0.72}$ M_⊙ yr⁻¹.
• 3.  
  Vinf: wind terminal speed V_∞ = 250, 500, and $1000\times {\hat{Z}}^{0.13}$ km s⁻¹.
• 4.  
  Rsc: cluster radius R_sc = 1 pc.
• 5.  
  age: current age t = 1 Myr.
• 6.  
  Mstar: total stellar mass M_⋆ = 2.05 × 10⁶ M_⊙.
• 7.  
  logLion: ionizing luminosity $\mathrm{log}{L}_{\mathrm{ion}}$ (erg s⁻¹).
• 8.  
  Namb: ambient density n_amb = 1, 10, 10², and 10³ cm⁻³.
• 9.  
  Tamb: mean ambient temperature T_amb determined by CLOUDY.
• 10.  
  Rmax: maximum radius ${R}_{\max }$ (pc) for the surface brightness integration.
• 11.  
  Raper: aperture radius R_aper (pc) for the total luminosity integration.
• 12.  
  Rshell: shell outer radius R_shell (pc).
• 13.  
  Rstr: Strömgren radius ${R}_{\mathrm{str}}$ (pc) determined by CLOUDY.
• 14.  
  Rbin: bubble inner radius R_b,in (pc).
• 15.  
  Rbout: bubble outer radius R_b,out (pc) or shell inner radius.
• 16.  
  Tbubble: median temperature T_bubble of the hot bubble (region (b) in Figure 2).
• 17.  
  Tadi: median temperature T_adi,med of the expanding wind predicted by the adiabatic solutions (region (a) in Figure 2).
• 18.  
  Twind: median temperature T_w,med of the expanding wind calculated by MAIHEM with the radiative solutions (region (a) in Figure 2).
• 19.  
  logUsp: dimensionless ionization parameter $\mathrm{log}{U}_{\mathrm{sph}}$ in a spherical geometry calculated by CLOUDY, defined as ${U}_{\mathrm{sph}}\equiv Q({{\rm{H}}}^{0})/4\pi {R}_{\mathrm{str}}^{2}{n}_{{\rm{H}}}c$, where Q(H⁰) is the total number of hydrogen-ionizing photons emitted per second, n_H the hydrogen density, and c the speed of light.
• 20.  
  thin: optically thin (1) or thick (0) model.
• 21.  
  mode: the cooling/heating radiative/adiabatic modes: 1 (AW: adiabatic wind), 2 (AB: adiabatic bubble), 3 (AP: adiabatic, pressure-confined), 4 (CC: catastrophic cooling), 5 (CB: catastrophic cooling bubble), and 6 (CP: catastrophic cooling, pressure-confined).
• 22.  
  H_1_1216, H_1_6563, ..., Ar_5_7006: integrated luminosities of the emission lines Lyα λ1216, Hα λ6563, ..., [Ar v] λ7006, respectively.