The Ionization and Dynamics of the Makani Galactic Wind

The detection of multiple emission lines throughout the Makani nebula allows us to constrain the spatially-resolved excitation of the wind using standard line ratios (Figure 5; Baldwin et al. 1981; Veilleux & Osterbrock 1987; Cid Fernandes et al. 2010).

Beginning with the nucleus, ap1, we that find that the broad component dominates its line flux and is securely in the AGN region of these diagrams. This is consistent with previous results (Perrotta et al., 2021). We extend this conclusion to the [O I]/H$\alpha$ vs. [O III]/H$\beta$ and [O I]/H$\alpha$ vs. [O III]/[O II] diagrams, based on our detection of [O I]. Perrotta et al. (2021) interpreted these line ratios as ionization from fast shocks; we discuss this in detail in Section 4. Makani lacks of evidence for an AGN in the rest-frame UV, X-ray, radio, or infrared (Diamond-Stanic et al., 2012; Sell et al., 2014; Petter et al., 2020; Whalen et al., 2022). The rest-frame optical-UV continuum of Makani shows no indication of non-stellar emission, and a stellar population fit to the SED is perfectly consistent with the data (Rupke et al., 2019). Though [Ne V] emission is present, its luminosity is low compared to AGN (Rupke et al., 2019). We thus discount AGN photoionization as the origin of these line ratios.

The narrow component of ap1 is instead likely ionized by stars, with some possible contribution from shocks based on its location in the composite region of the [N II]/H$\alpha$ vs. [O III]/H$\beta$ diagram.

The two apertures with $r\sim 10$ kpc (ap2 and ap3) combine information from two velocity components, one broad and one narrow. Ap2 has line ratios similar to those from total fluxes in ap1, with higher low-ionization line ratios ([O I]/H$\alpha$) and lower ionization parameter [O III]/[O II]. Apertures ap3, ap4, and ap5 have similar [O I]/H$\alpha$ and [O III]/[O II] as ap2, but with lower high-ionization line ratios ([O III]/H$\beta$), larger [S II]/H$\alpha$, and slightly smaller [N II]/H$\alpha$. Two other apertures with significant detections across several lines (ap6 and ap8) lie in roughly the same locations as the other apertures, but these apertures have low S/N in the weakest lines so their positioning is more uncertain. The lines other than [O II] and H$\alpha$ in ap7 and ap9 are not strong enough to accurately place them in these diagrams, though ap9 does appear in several with large error bars.

Taken as a whole, the apertures other than ap1 lie largely in the composite region of the [N II]/H$\alpha$ diagram and in the low-ionization nuclear emission-line region (LINER) area in the other three diagrams. An exception is aperture ap2, which lies in between the properties of ap1 and ap3. This aperture is in the direction of the most highly blueshifted molecular and ionized gas (Rupke et al., 2019). It must therefore include more flux from this higher-ionization wind component (mixed in with a lower-ionization component) than ap3, which is at the same radius as ap2.

We detect the weak emission lines [O III] 4363 Å and [N II] 5755 Å in ap1, and [O III] 4363 Å in ap2. The [O III] 4363/5007 line flux ratios in ap1 and ap1.2 are -1.1 dex. Ratios this high are observed in LINERs, implying high temperatures ($T\gtrsim 2\times 10^{4}$ K) (Nagao et al., 2001; Molina et al., 2018). Similarly, the ratio [N II] 5755/6548 is in the range -1.3 to -1.4, also consistent with $T\sim 2\times 10^{4}$ K.

Makani may be the highest-redshift system in which interstellar Na I D absorption and/or emission has been seen in a galactic wind. Resonant Na I D absorption is detected in the host galaxy (ap1) with a low equivalent width (rest-frame $W_{\mathrm{eq}}$ of 0.37${}^{+0.07}_{-0.05}$ Å), but blueshifted by -342 km s^-1 (Figure 3). The outflow was not previously detected in resonant absorption in Mg II 2796, 2803 Å or Mg I transitions in the near-UV (Rupke et al., 2019). However, there is residual, weak Fe II 2585 Å absorption (Rupke et al., 2019; Perrotta et al., 2021).

Previously, redshifted resonant emission in Mg II, Mg I, and Fe II^∗ was detected in Makani (Rupke et al., 2019). The Mg II is spatially resolved into an $r\sim 20$ kpc nebula. We now detect spatially-extended Na I D emission in ap1 and ap2, with fluxes and rest-frame equivalent widths of $(3.2\pm 0.2)\times 10^{-17}$ erg s^-1 cm^-2 and $W_{\mathrm{eq}}=-0.49\pm 0.05$ Å in ap1 and $(2.0\pm 0.2)\times 10^{-17}$ erg s^-1 cm^-2 and $W_{\mathrm{eq}}=-4.4\pm 0.5$ Å in ap2.

Several lines of evidence suggest that resonant absorption in ap1 is filling in the emission, reducing the observed $|W_{\mathrm{eq}}|$ and increasing $|v-v_{\mathrm{sys}}|$ of both absorption and emission. First, a fit to Fe II 2585 Å gives a velocity at maximum depth of -100 km s^-1 and a covering factor of 0.12, assuming it is optically thick. Thus $v_{\mathrm{NaI}}-v_{\mathrm{FeII}}\sim-200$ km s^-1. Furthermore, if the Na I D and Fe II absorption lines are optically thick, the ratio of the observed covering factors is $C_{\mathrm{f,FeII}}/C_{\mathrm{f,NaI}}\sim 3-4$, suggesting that Na I D is too shallow.

Second, neutral gas absorption in outflows is closely connected with the foreground dust column (Rupke & Veilleux, 2015; Rupke et al., 2021). Though there is significant scatter, the observed $E(B-V)=0.6-1.0$ in ap1 and ap2 (Figure 4) are typically associated with 5–10$\times$ higher $W_{\mathrm{eq}}$ on average.

Finally, a comparison to the radial surface-brightness profile of Mg II (Figure 6) shows that the Na I D surface brightness declines much more slowly than that of Mg II. If this is due to infilling of emission by foreground absorption, Na I D should intrinsically be brighter in ap1 by a factor of a few ($0.6\pm 0.2$ dex). This in turn would raise the absorption equivalent width by the same factor and the covering factor by some (smaller) amount. This reduction of Na I D flux due to infilling is observed in other spatially-resolved observations (Rupke & Veilleux, 2015).

Due to this infilling, we do not attempt to estimate the column density, and thus the mass and dynamics, of the neutral phase of the wind.

Makani contains a compact ($r\sim 400$ pc) stellar core (Sell et al., 2014; Rupke et al., 2019; Diamond-Stanic et al., 2021) with a young stellar population (Rupke et al., 2019). The low-velocity component of the innermost ($r\lesssim 5.5$ kpc) aperture, ap1.1, thus contains most or all of the star formation activity. Based on the observed line ratios (Section 3.1 and Perrotta et al. 2021), this component is photoionized by the young stellar population. The H$\alpha$ luminosity of this component is $2.1\times 10^{42}$ erg s^-1 (Table 2), which is 40% of the total luminosity in ap1 (i.e., the host galaxy). We note that this luminosity, however, is larger than the difference between the total nebular luminosity (Section 3) and wind luminosity (Section 4.4). This difference is $L_{\mathrm{H\alpha}}^{\mathrm{starburst}}=L_{\mathrm{H\alpha}}^{\mathrm{tot}}-L_{\mathrm{H\alpha}}^{\mathrm{wind}}=(9.4-8.0)\times 10^{42}$ erg s${}^{-1}=1.4\times 10^{42}$ erg s^-1, 30% lower than the value measured from ap1.1. We take this difference as an indication of the uncertainties in our methodology, which includes assigning fluxes among the starburst and wind Episodes I and II, bootstrapping from the KCWI [O II] data from the ESI meaurements, and the absolute flux calibration of the current data (Section 2). In what follows we assume the average of these two values, $L_{\mathrm{H\alpha}}^{\mathrm{starburst}}=1.8\times 10^{42}$ erg s^-1.

Makani has a star formation rate (SFR) of 224–300 M_⊙ yr^-1, as inferred from radio and infrared data (Petter et al., 2020). The radio emission is concentrated in a 1″ point source coinciding with ap1. Using a standard calibration (Moustakas et al., 2006), we compute from $L_{\mathrm{H\alpha}}^{\mathrm{starburst}}$ a SFR of 14 M_⊙ yr^-1. This is a factor of 16 smaller than the lowest previous estimate from IR data. A similar discrepancy is seen in other galaxies in the parent galaxy sample of Makani (Edmonds et al., 2022). It may in fact be evidence of the stellar feedback in Makani shutting off star formation, as H$\alpha$ traces star formation over the shortest timescales $\lesssim$6 Myr (Calzetti, 2013). This is just below the estimated 7 Myr estimate of the most recent burst of star formation. An earlier burst occurred 400 Myr ago (Rupke et al., 2019), which is unlikely to be reflected in current SFR estimates. The observed discrepancy could also result in part from LyC leakage (Moustakas et al., 2006), to which compact starbursts with strong feedback like Makani may be susceptible (Perrotta et al., 2021).

The bright, extended Mg II emission in Makani is consistent with some LyC escape (Chisholm et al., 2020; Xu et al., 2022). The observed Mg II flux ratio in Makani, based on KCWI spectra, is $R\equiv F_{2796}/F_{2803}=1.00\pm 0.08$. Using the optically-thin scenario, or equivalently the clumpy scenario with optically-thin channels, from Chisholm et al. (2020) yields an escape fraction in the Mg II 2796 Å line of 25%. Using the line fluxes from ap1.1 in the current data and the method of Xu et al. (2022), we can make a second, independent measure of $f_{\mathrm{esc}}(2796)$. We apply the highest-metallicity model relating the intrinsic values of [O III]5007/[O II] and Mg II2796/[O III]5007 (Xu et al., 2022) to yield $F_{2796}(\mathrm{intrinsic})=2.92\times 10^{-16}$ erg s^-1 cm^-2 in the optically-thin, density-bounded case. In a 1$\farcs$5-diameter nuclear KCWI aperture, roughly corresponding to the footprint of ap1.1, we measure $F_{2796}(\mathrm{observed})=0.70\times 10^{-16}$ erg s^-1 cm^-2. Then $f_{\mathrm{esc}}(2796)=F_{2796}(\mathrm{obs})/F_{2796}(\mathrm{int})=24\%$, identical to the first estimate. Comparing to observations of other Mg II emitters, this value of $f_{\mathrm{esc}}(2796)$ suggests that $f_{\mathrm{esc}}(\mathrm{Lyc})$ could be a few percent or higher.

Finally, it may also be the case that the bulk of the star formation is obscured by larger columns than traced by optical light, as seen in nearby compact mergers (Spoon et al., 2007). Future rest-frame far-ultraviolet and mid-infrared spectra of Makani and its parent sample could distinguish among these possibilities.

The oxygen abundance of this component is 12$+$log(O/H) $\sim$ 9.1 (or $\sim$2$\times$ Solar; Asplund et al. 2021), based on standard [N II]/[O II] and $R_{23}$ strong-line calibrations (Kewley & Dopita, 2002). For component 1 of ap1, [N II]/[O II] $=-0.12$ and $R_{23}=-0.13$. This is consistent with the high mass of the galaxy (Rupke et al., 2019; Diamond-Stanic et al., 2021). The actual abundance could be closer to Solar, as strong-line calibrations may overestimate the true gas abundance by a factor of $\sim$2 (Kewley & Ellison, 2008). The weak lines [O III] 4363 Å and [N II] 5755 Å could in principle also provide a weak-line measurement based on electron temperature. However, the inferred temperatures are much too high for high-metallicity H II regions, suggesting that component blending or simply an overestimate of the ap1.1 component in these weak lines is producing these high ratios.

Line ratios are commonly used to distinguish between stellar and shock ionization in galactic winds (Lehnert & Heckman, 1996; Veilleux & Rupke, 2002; Sharp & Bland-Hawthorn, 2010). Aside from ap1.1, the line ratios observed in Makani (Figure 5) are inconsistent with those predicted by stellar photoionization models (Kewley et al., 2001; Byler et al., 2017). An obvious place to turn instead is shock models, both due to their typical consistency with LINER-like line ratios (e.g., Dopita et al., 1997, among numerous examples) and because of the extreme outflows we detect in the Makani nebula.

We first compare the position of apertures ap3, ap4, and ap5 in Figure 5 to the low-velocity range of fast shock models (Allen et al., 2008). We assume solar metallicity, and the benchmark model pre-shock density $n_{\mathrm{H}}=1$ cm^-3 (their model series M_n1). In the [N II]/H$\alpha$ diagram, shock velocities $v_{\mathrm{s}}\sim 200$ km s^-1 and low magnetic field strengths $B/n^{1/2}<0.5$ $\mu$G cm^-3/2 can reproduce the results, whether the emission comes from the shock only or from the shock and precursor. The precursor is the gas that is about to be impacted by the shock, and which in this case is fully pre-ionized by the shock radiation. For [S II]/H$\alpha$, the same low $v_{s}$ and low $B$ fit best in the shock-only grid, but the shock$+$precursor grid does not overlap the observed values. The [O I]/H$\alpha$ grids allow for shock or shock$+$precursor solutions, though with slightly higher velocity, $v_{\mathrm{s}}\sim$250 km s^-1.

The fast shock models cannot match the low observed ionization parameters from [O III]/[O II], even in the lower-ionization shocked region. This may point to only partial pre-ionization of the precursor, as assumed in slow shock models (Rich et al., 2010; Dopita & Sutherland, 2017), or additional post-shock physics. The slow-shock models are consistent with the observed [S II]/H$\alpha$ and [O I]/H$\alpha$ at the highest model velocities, $v_{\mathrm{s}}\sim 200$ km s^-1, but produce values of [N II]/H$\alpha$ that are too high at the same velocities.

The fast shock models are mostly consistent with the observed values of [O II]/H$\alpha=1-5$ (Figure 4). For shocked gas only, the model ratio is 2–4 for $v_{\mathrm{s}}=100-300$ km s^-1, with a minimum around 200 km s^-1. Including precursor, the model ratio is slightly smaller, in the 1.8–2.8 range. The observed values are much larger than those in nearby star-forming galaxies, for which the maximum observed [O II]/H$\alpha$ is about 2 (Moustakas et al., 2006). It is also larger than observed in the diffuse ionized gas in extraplanar regions of nearby galaxies (Jones et al., 2017).

We estimate the emitted H$\beta$ fluxes, $F_{\mathrm{H\beta}}^{\mathrm{shock}}$, across the shock in the context of these models by dividing the luminosity in each aperture (Table 2) by the projected physical area of each aperture in cm^-2 (Table 1). This assumes the shock is plane-parallel to the line-of-sight, and so may be an underestimate if the aperture is tracing the edge of the wind. However, given the large size of the wind, projection effects are probably not significant for the inner apertures (ap1–ap5). The resulting shock fluxes $F_{\mathrm{H\beta}}^{\mathrm{shock}}$ in ap3–ap5 are approximately 10^-4.0 to 10^-3.7 erg s^-1 cm^-2. We then compare these observed fluxes to those predicted by the the $Z=Z_{\odot}$ models M_n1 and V_n10 from Allen et al. (2008). We use the shock velocity and magnetic field parameter estimated from comparing the observed line ratios to the same models, $v_{\mathrm{s}}=200$ km s^-1 and $B/n^{1/2}<0.5$ $\mu$G cm^-3/2. (The result is not sensitive to $B$.) Finally, interpolating between the model grids, the extrapolated pre-shock density from fast models is then $n_{\mathrm{H}}=2-3$ cm^-3, since the flux scales linearly with $n_{\mathrm{H}}$ (Dopita & Sutherland, 1995). At larger radius, the H$\beta$ fluxes decrease to the range $10^{-4.3}$ to $10^{-5.7}$ erg s^-1 cm^-1, which suggests an increase in projection effects that mean these are lower limits; or a lower velocity and/or density at these radii. Lower velocity would be consistent, for instance, with the lower line ratios in ap6 and ap8 compared to ap3–ap5, but it seems probable that the gas density would also decrease with radius.

At $r=$ 0–5 kpc, the high-velocity outflow component ap1.2 shows elevated [O III]/H$\beta$, [N II]/H$\alpha$ and [O III]/[O II]; and lower [O I]/H$\alpha$ and [S II]/H$\alpha$. The H$\beta$ flux is also higher at 10^-2.65 erg s^-1 cm^-1. The line ratios place this broad component in the AGN regions of the excitation diagrams. However, $n_{\mathrm{H}}\sim 10$ cm^-3 shock$+$precursor models from Allen et al. (2008) (their models V_n10) are entirely consistent with the observed line ratios if $v_{\mathrm{s}}$ is in the range 300–400 km s^-1 and $B/n^{1/2}<0.5$ $\mu$G cm^-3/2. These models also reproduce the observed ratios [Ne V] 3426/H$\beta=(-1.33\pm 0.20)$ dex and [Ne V] 3426/[Ne III] 3869 $=(-0.45\pm 0.22$) dex. The higher ionization state is due to the harder radiation field produced by the higher-velocity shocks. This harder radiation also produces high temperatures in the precursor and post-shock regions. The observed [O III] 4363/5007 ratio in ap1.2 is -1.1 dex, consistent with these high temperature ($T\gtrsim 2\times 10^{4}$ K; Section 3). Such ratios are commonly observed in LINERs (Nagao et al., 2001; Molina et al., 2018). However, the models cannot perfectly reproduce the combined [O III] 4363/5007 and [O III]/H$\beta$ ratios, an ongoing problem (Dopita & Sutherland, 1995; Allen et al., 2008).

In summary, the fast shock models of Allen et al. (2008) are very consistent with the observed strong line ratios and fluxes throughout the nebula. In ap1.2 (the outflow at radii 0–5 kpc), high shock velocities (300–400 km s^-1) and an ionized precursor are needed. These shock velocities align with the observed $\langle\sigma\rangle_{\mathrm{II}}=400$ km s^-1 (Rupke et al., 2019). Farther from the host, the required shock velocities are lower, 200 km s^-1, but these still align with the observed linewidths $\langle\sigma\rangle_{\mathrm{I}}=200$ km s^-1. The precursor emission contributes fractionally less at these velocities. In the next section, we apply the shock models to constrain the wind density.

We can combine the observed line flux ratios [S II]6716/6731 and [O II]3729/3726 with the Allen et al. (2008) shock models to constrain the ambient, pre-shock and post-shock densities in the wind and extended nebula. The high-velocity Episode II wind seen in ap1–ap3 has a high electron density $n_{e}=3000$ cm^-3 (1$\sigma$ range 1200–10,000 km s^-1) from the [S II] ratio of ap1.2. This density is also consistent with the [O II] ratio. The shock model that best fits the observed line ratios and shock fluxes in ap1.2–ap3 has a pre-shock density of $n_{\mathrm{H}}\sim 10$ cm^-1, as we discuss above in Section 4.2. However, compression in the shock raises the density in the post-shock region up to a factor of $n_{\mathrm{post}}/n_{\mathrm{pre}}=(8\pi\mu m_{\mathrm{H}})^{1/2}v_{\mathrm{shock}}/(B/n^{1/2})$ (Dopita & Sutherland, 1996). For the estimated velocities of 300–400 km s^-1 and low magnetic field parameter $0.5~{}\mu$G cm^-3/2, the maximum value of log($n_{\mathrm{post}}/n_{\mathrm{pre}})$ is 2.75. This translates to a post-shock density of at most $n_{\mathrm{post}}\sim 6000$ cm^-3 for $n_{\mathrm{pre}}=10$ cm^-3. (The compression can also be seen in Figures 3–6 of Allen et al. (2008).)

The correspondence between the observed density and predicted post-shock density suggests that the observed [S II] and [O II] emission is arising primarily in the compressed post-shock region. The line ratios, however, appear to require an ionized precursor. That the post-shock region seems to dominate the flux of these lines could point to additional physics occurring behind the shock front, perhaps due to the accumulation of swept-up interstellar gas or to the driving mechanism of the wind.

At larger radii, in the extended, Episode I nebula, the line ratios and fluxes point to much lower ionized gas densities $n_{e}\lesssim 10$ cm^-3. The shock models are consistent with low pre-shock densities of $n_{\mathrm{H}}\lesssim 2-3$ cm^-3, with this value possibly decreasing with increasing radius. The corresponding maximum post-shock density is $n_{\mathrm{H}}\lesssim 600$ cm^-3. In contrast to the high-velocity, high-density Episode II gas, this implies that the [S II] ratio is primarily tracing the lower-density precursor gas. The line ratios are, for the most part, consistent with either shock-only or shock$+$precursor models.

The 200 and 400 km s^-1 models with $n_{\mathrm{H}}=1-10$ cm^-3 (models M_n1 and V_n10) predict that just a little over half of the H$\beta$ luminosity is produced in the post-shock gas, while the other fraction arises in the precursor, with little dependence on magnetic parameter $B/n^{1/2}$ (Allen et al., 2008). Thus, in the calculations that follow, we assume complete ionization in each shock region and divide the flux evenly between these low- and high-density regions. Using the estimated pre-shock densities and calculated maximum compression as a guideline, we assume $n_{e}=10$ and 1000 cm^-3 pre- and post-shock for Episode I, and $n_{e}=1$ and 100 cm^-3 for Episode II. The post-shock densities could be even higher, but higher values do not significantly change the results, as the lower-density pre-shock gas dominates the mass.

As the KCWI observations did not cover strong H recombination lines, we were previously unable to measure the mass of the ionized Makani nebula. Here we use our bootstrapped estimate of the spatially-resolved H$\alpha$ luminosity of the wind (Section 3) to estimate the mass and dynamics of both the Episode I and II winds.

We start with the Voronoi-binned, KCWI [O II] map and spatially delineate the Episode I and II winds using a velocity cut at $v_{98\%}=-700$ km s^-1 (Rupke et al., 2019). To do this, we characterize multi-component line fits in the KCWI data with the cumulative velocity distribution. We calculate velocities that enclose a particular percentage p of flux integrated from the red side of the line, $v_{\mathrm{p\%}}$. (E.g., $v_{98\%}$ is the velocity at which 98% of the line flux is redshifted from this velocity, and so characterizes the maximum blueshift observed.) We compute the H$\alpha$ flux and luminosity in each bin, applying our model for ([O II], observed)/(H$\alpha$, intrinsic) vs. radius and using the mean radius in each Voronoi bin. From the inner $r<5.5$ kpc, we also conservatively remove 40% of the flux in each bin for the contribution from star-forming gas (Section 4.1).

As discussed in Section 4.3, we use the observed line ratios and shock models to infer how much of the observed luminosity arises from regions of different density. Using these luminosities and densities, we then calculate the ionized gas mass in each bin. This ionized density model differs from many estimates in the literature, which typically use only the density computed from line ratios. However, the shock models provide further information on the density structure that could yield a more accurate estimate. Given that $M^{\mathrm{HII}}\sim n_{e}^{-1}$, the masses are largely dependent on the lower, pre-shock density in our model. This lower density is mostly consistent with the measured densities (Section 3). However, the pre-shock density in the Episode II shock model is lower than the density measured in the densest part of the wind, in ap1.2, by a factor of $\sim$100. Assuming a constant density wind with the higher density measured from [S II] and [O II] in ap1.2 would result in a mass outflow rate that is lower by a similar factor.

The resulting H$\alpha$ luminosities and gas masses are listed in Table 3. The Episode I gas is fainter but contains more mass because of our density model. Together, the total ionized gas mass of $6.4^{+5.2}_{-2.6}\times 10^{9}$ M_⊙ is comparable to the integrated CO mass of $1\times 10^{10}$ M_⊙ (Rupke et al., 2019). The detected CO is much less extended than the ionized gas (reaching a radius of only 20 kpc; Figure 6). The Episode II ionized gas mass is also comparable to the CO mass of $(2.4\pm 0.6)\times 10^{9}$ M_⊙ in a similar velocity range (Rupke et al., 2019).

To compute outflow rates (Table 3), we use two methods. First, we assume a single-radius wind for each episode, using $R_{I}=40$ kpc and $R_{\mathrm{II}}=15$ kpc (Method 1). We use both the central velocity $v_{50\%}$ (Method 1a) and maximum velocity $v_{98\%}$ (Method 1b) to compute the wind properties. We deproject the velocity (either $v_{50\%}$ or $v_{98\%}$) in each bin to compute the 3D $v_{\mathrm{rad}}$, using the projected galactocentric distance of the bin $R_{\mathrm{obs}}$ and assuming that it lies on the wind surface at the assumed 3D radius $R_{\mathrm{wind}}$: $v_{\mathrm{rad}}=v_{obs}/\mathrm{cos}[\mathrm{sin}^{-1}(R_{\mathrm{obs}}/R_{\mathrm{wind}})]$. We then compute the ballistic flow time for each bin, $t_{\mathrm{flow}}=v_{\mathrm{rad}}/R_{\mathrm{wind}}$ and divide the mass in that bin by $t_{\mathrm{flow}}$. We follow a similar procedure for momentum and energy, adding in velocity dispersion for energy (following Rupke et al. 2005).

As an alternative method we assume the stellar population ages $t_{*}$ for Episodes I and II, 400 and 7 Myr, represent the times when all of the gas in each episode left the starburst (Method 2). We then divide the ionized mass by these times for each episode to compute $dM/dt=M/t_{*}$. This method is effectively a time-averaged outflow rate, where the time average is over the lifetime of the outflow. We use it as a sanity check on Method I.

For Episode I, Method 1a (using $v_{50\%}$) yields a result consistent with Method 2, as the mean deprojected radial velocity $\langle v_{\mathrm{rad}}\rangle=90$ km/s is close to the velocity inferred from $R_{\mathrm{I}}/t_{\mathrm{*,I}}=40~{}\mathrm{kpc}/400~{}\mathrm{Myr}=97$ km/s. In other words, most of the gas is at a velocity consistent with the outer radius of the nebula and the ballistic flow velocity, if this gas left the starburst 400 Myr ago. The gas that began at higher velocity has either slowed or since escaped to much larger radius. The remaining detected amounts of high-velocity gas do not contribute significantly to the mass, or mass outflow rate, of the outflow. Using $v_{98\%}$ to trace the radial velocity of the bulk of the gas thus overestimates the flow rate.

For Episode II, Method 1b (using $v_{98\%}$) is consistent with Method 2, due again to comparable $\langle v_{\mathrm{rad}}\rangle=2279$ km/s and $R_{\mathrm{II}}/t_{*,\mathrm{II}}=15~{}\mathrm{kpc}/7~{}\mathrm{Myr}=2090$ km/s. Method 1a, in contrast, yields a flow rate 15–20$\times$ lower. In this case, while the fastest gas has flowed to the largest observed radius, the lower velocity gas has not had time to do so. Thus, the assumption of a single radius for the inner wind is probably incorrect, and a flow rate computed from the flow time in each spaxel, $t_{*}$ (i.e., Method 2), is more correct. This intepretation implies that the inferred radius, $R_{\mathrm{II}}=v_{\mathrm{rad}}/t_{*,\mathrm{II}}$, is velocity-dependent. Lower-velocity spaxels will have smaller radius, and our single-radius deprojection will underestimate the radial velocity. For both of these reasons, Method 1a produces a significant underestimate for the Episode II wind.

We conclude that Methods 1a and 2 are closer to the correct flow rates for Episode I, while Methods 1b and 2 are closer to the correct flow rates for Episode II. The inferred mass outflow rates are thus in the range 10–13 M_⊙ yr^-1 for Episode I and 151–186 M_⊙ yr^-1 for Episode II. Again, $dM/dt_{\mathrm{II}}$ for the ionized gas is comparable to the CO value, 245 M_⊙ yr^-1. In most nearby compact starbursts, $dM/dt$ in the ionized gas is much smaller than in the molecular gas, though with significant scatter (e.g., Rupke & Veilleux, 2013; Fluetsch et al., 2021). However, the canonical density model in these calculations is inferred from optical line ratios; we combine the line ratio densities with the shock model density structure, which raises the inferred masses. In Perrotta et al. (2022, in prep.), we will show that measurements of $dM/dt$ from Mg II and Fe II absorption lines in other compact starbursts from the Tremonti et al. (2007) sample are consistent with mass outflow rates of order $10^{2-3}$ M_⊙ yr^-1.

Note that we choose the single wind radius for each episode based on two factors. The first is the observed wind morphology (Rupke et al., 2019). We refine the radius to approximately match the average deprojected velocity $\langle v_{\mathrm{rad}}\rangle$ and the velocity $R/t_{*}$ for the two methods (1a or 1b and 2) that produce the best agreement in $dM/dt$ for each episode. The resulting $dM/dt$ is not sensitive to a 30% change in $R_{\mathrm{wind}}$.

We have shown that the line emission from the wind is most consistent with shock models (Section 4.2). We nonetheless consider if the radiative energy from the starburst is capable of powering the H$\alpha$ luminosity of the nebula.

We first assume the current SFR as given by IR/radio tracers, rather than H$\alpha$. A $\lesssim$6.5 Myr instantaneous burst or a $\gtrsim$7 Myr continuous burst with star formation rate of 224–300 $M_{\odot}$ yr^-1 produces $Q(\mathrm{H}^{0})=(1.9-2.6)\times 10^{55}$ ionizing photons s^-1 (Leitherer et al., 1999). The ratio of H$\alpha$ photons to ionizing photons is $\alpha_{\mathrm{H\alpha}}^{\mathrm{eff}}/\alpha$, or the ratio of the effective to total recombination coefficients. For Case B at $T=10^{4}$ K and $n_{e}=100$ cm^-2, $\alpha_{\mathrm{H\alpha}}^{\mathrm{eff}}/\alpha_{\mathrm{B}}=\alpha_{\mathrm{H\beta}}^{\mathrm{eff}}/\alpha_{\mathrm{B}}\times(j_{\lambda}\lambda)_{\mathrm{H\alpha}}/(j_{\lambda}\lambda)_{\mathrm{H\beta}}=0.452$, where $j_{\lambda}$ are the line emissivities and we use the emissivity and recombination coefficient tabulations of Hummer & Storey (1987). This yields an H$\alpha$ luminosity predicted from the radio/IR SFR of $L_{\mathrm{H\alpha}}^{\mathrm{SFR}}=(2.6-3.5)\times 10^{43}$ erg s^-1, which is $3-4$$\times$ larger than the observed luminosity of the nebula, $L_{\mathrm{H\alpha}}^{\mathrm{tot}}(\mathrm{H}\alpha)=9.4\times 10^{42}$ erg s^-1 (Section 3).

Photoionizing the entire nebula would require a significant fraction of the Lyman continuum to leak outside of the inner starburst region. We can infer from the predicted H$\alpha$ luminosity $L_{\mathrm{H\alpha}}^{\mathrm{SFR}}$ and the H$\alpha$ emission attributed to the starburst, $L_{\mathrm{H\alpha}}^{\mathrm{starburst}}=1.8\times 10^{42}$ erg s^-1 (Section 4.1), that a fraction $L_{\mathrm{H\alpha}}^{\mathrm{starburst}}/L_{\mathrm{H\alpha}}^{\mathrm{SFR}}\sim 0.05-0.07$ of the LyC is required to ionize the gas consistent with stellar photoionization (ap1.1). To energize the H$\alpha$ emisison in the wind ($L_{\mathrm{H\alpha}}^{\mathrm{wind}}=8.0\times 10^{42}$ erg s^-1; Table 3) would require an additional 23–31% of the LyC emitted by the starburst. The total escape fractions for the inner starburst and entire nebula, respectively, would then be of order 90% and 65%, which are uncomfortably high for a galaxy of this mass and ionization parameter (Izotov et al. 2022). Thus, if the nebula were photoionized by the starburst rather than shocks, significant absorption of LyC by dust or a very recent SFR that is much lower would be necessary (Section 4.1).

The shock models also make energetic predictions about the total flux radiated in the shocks to which we can compare the data. The mechanical energy flux that powers the shock is completely radiated by the cooling gas (Dopita & Sutherland, 1996). Assuming H$\alpha$/H$\beta=2.86$, the ratio of the total radiated luminosity $L_{\mathrm{rad}}$ to the observed H$\alpha$ luminosity is given by (Dopita & Sutherland, 1995)

Here we compare the energy produced in the shocks to the mechanical energy observed in the warm, ionized and molecular phases of the wind episodes. We also look at how the mass, momentum and energy in the wind could arise from the starburst itself.

Observations of molecular gas and Mg II 2796, 2803 Å emission to radii 20–25 kpc point to the presence of significant H I and H₂ in the inner, Episode II wind. The presence of these two tracers also indicates significant amounts of dust out to similar radii. However, these dusty gas phases appear only at the inner edge of the Episode I wind (Rupke et al., 2019).

Two new lines of evidence confirm this picture. First, Na I D is detected in two apertures out to $r=10-15$ kpc (Figure 6). The presence of Na I D has long been closely associated with dust extinction (see, e.g., Rupke et al., 2021, and references therein), due to dust being required to keep Na in the neutral phase. Second, we detect substantial gas extinction, $E(B-V)\sim 0.5$, out to $r=20-25$ kpc (Figure 4).

Such dusty, neutral outflows on $\sim$10 kpc scales are commonly observed in compact starbursts in nearby mergers (Rupke & Veilleux, 2013). Frequently, high spatial resolution observations show filamentary dust structures associated with these outflows. Future high-resolution observations of Makani with the James Webb Space Telescope may resolve these dust structures, as well as more directly determine the radial extent of dust in the wind.