Non-equilibrium Ionization States Within Galactic Outflows: Explaining Their O VI and N V Column Densities

TL;DR

This study uses hydrodynamic simulations with non-equilibrium chemistry to explain observed ionization states in galactic outflows, revealing departures from equilibrium assumptions and matching observed ion ratios.

Contribution

It introduces a comprehensive non-equilibrium ionization model for galactic outflows, challenging the common equilibrium assumption in observational analysis.

Findings

Ionization states often depart from equilibrium predictions.

Models reproduce observed N V to O VI ratios.

Simulations match observed O VI column densities.

Abstract

We present a suite of one-dimensional spherically-symmetric hydrodynamic simulations that study the atomic ionization structure of galactic outflows. We track the ionization state of the outflowing gas with a non-equilibrium atomic chemistry network that includes photoionization, photo-heating, and ion-by-ion cooling. Each simulation describes a steady-state outflow that is defined by its mass and energy input rates, sonic radius, metallicity, and UV flux from both the host galaxy and meta-galactic background. We find that for a large range of parameter choices, the ionization state of the material departs strongly from what it would be in photo-ionization equilibrium, in conflict with what is commonly assumed in the analysis of observations. In addition, nearly all the models reproduce the low N V to O VI column density ratios and the relatively high O VI column densities that are observed.

Tables

Name	$\dot{M}$	$\dot{E}$	$\alpha$	$\beta$	Voutflow	n	T	Z	UV	SFR
	[M⊙/yr]	[1050erg/yr]			[km/s]	[cm-3]	[106 K]	[Z⊙]		[M⊙]
Nominal	10	1.0	1.0	1.0	500	4.6	18	1.0	✓	10
UV0	10	1.0	1.0	1.0	500	4.6	18	1.0		10
Z05	10	1.0	1.0	1.0	500	4.6	18	0.5	✓	10
Z20	10	1.0	1.0	1.0	500	4.6	18	2.0	✓	10
M05	5	1.0	1.0	0.5	707	1.7	36	1.0	✓	10
M02	2	1.0	1.0	0.2	1100	0.4	91	1.0	✓	10
M05E05	5	0.5	0.5	0.5	500	2.3	18	1.0	✓	10
M50E50	50	5.0	0.5	0.5	500	23	18	1.0	✓	100
SLow	1.5	0.47	1.0	0.3	885	0.4	57	1.0	✓	5
SHigh	12	1.7	0.9	0.6	594	4.7	26	1.0	✓	20

Equations

$$\frac{1}{r^2}\frac{d}{dr}(\rho u{r}^2)=q_m,$$

$$\rho u\frac{du}{dr}=-\frac{dP}{dr}-q_{m}u,$$

$$\frac{1}{r^2}\frac{d}{dr}\left[\rho u{r}^2\left(\frac{1}{2}{u}^2+\frac{\gamma }{\gamma -1}\frac{P}{\rho }\right)\right]=q_e,$$

$$q_m=\{\begin{array}{ll}\dot{M}/V& \text{if}r\le R_{\ast }\\ 0& \text{if}r>R_{\ast },\end{array}$$

$$\begin{gathered} q_e= \\ \{\begin{array}{ll}\dot{E}/V& \text{if}r\le R_{\ast }\\ -\Lambda n^2-\mathcal{H}n& \text{if}r>R_{\ast },\end{array} \end{gathered}$$

$$\Lambda n^2=\sum _i{\Lambda }_i{n}_i{n}_e,$$

$$\mathcal{H}n=\sum _i{\mathcal{H}}_i{n}_i,$$

$$\dot{E}=\dot{M}\left(1/2v_{\mathrm{outflow}}^2+3/2c_s^2\right),$$

$$c_s^2=v_{\mathrm{outflow}}^2=\dot{E}/2\dot{M}.$$

$$\Omega R_{\ast }^2\rho v_{\mathrm{outflow}}=\dot{M},$$

$$T_{\mathrm{outflow}}\approx v_{\mathrm{outflow}}^2{m}_H\bar{A}/k_B,$$

$$L_{\mathrm{outflow}}=L_{6,SB99}\left(\frac{{\dot{M}}_{SFR}\ast t_{\mathrm{evo}}}{10^6{M}_{\odot }}\right),$$

$$B_{\nu ,SB99}=L_{\mathrm{outflow}}\left(\frac{{\lambda }^2}{c}\right)\left(\frac{1}{16{\pi }^2{R}_{\ast }^2}\right),$$

$$B_{\nu }(R)=\left(B_{\nu ,SB99}+B_{\nu ,\mathrm{HM}}\right){\left(\frac{R_{\ast }}{R}\right)}^2,$$

Full text

galaxies: evolution – galaxies: starburst – methods: numerical – hydrodynamics – ISM: jets and outflows

††software: FLASH (Fryxell et al., 2000), yt (Turk et al., 2011), CLOUDY(Ferland et al., 2013)

1 Introduction

Galactic outflows are ubiquitous in intensely star-forming galaxies across all redshifts (Heckman et al., 1990; Lehnert & Heckman, 1996; Heckman et al., 2000; Pettini et al., 2002; Martin, 2005; Rupke et al., 2005; Veilleux et al., 2005; Weiner et al., 2009; Bordoloi, 2014; Bordoloi et al., 2016; Rubin et al., 2014; Heckman et al., 2015; Chisholm et al., 2016). These outflows are powered by stars through supernovae (e.g., Mac Low & Ferrara, 1999; Scannapieco et al., 2001; Mori et al., 2002; Scannapieco et al., 2002; Springel & Hernquist, 2003; Dalla Vecchia & Schaye, 2008; Creasey et al., 2013), stellar winds (e.g., Murray et al., 2011; Hopkins et al., 2012; Muratov et al., 2015; Hayward & Hopkins, 2017), and cosmic rays (e.g., Farber et al., 2018). They have dramatic effects on their host galaxies by decreasing metallicities (e.g., Tremonti et al., 2004; Oppenheimer et al., 2009; Davé et al., 2011; Lu et al., 2015; Agertz & Kravtsov, 2015), suppressing star-formation (e.g., Somerville & Primack, 1999; Cole et al., 2000; Scannapieco et al., 2001, 2002; Benson et al., 2003) and perhaps occasionally enhancing it (e.g., Scannapieco et al., 2004; Gray & Scannapieco, 2010, 2011a, 2011b; Bieri et al., 2016; Fragile et al., 2017; Mukherjee et al., 2018).

The nature of galactic outflows is very complex, and a complete picture is possible only when they are observed and modeled over the full range of temperatures – from X-ray observations of 107 to 108 K material (Martin, 1999; Strickland & Heckman, 2007, 2009), near-UV and optical observations of warm ionized gas at $\approx$104 (e.g., Pettini et al., 2001; Tremonti et al., 2007; Martin et al., 2012; Soto et al., 2012), and IR and submm observations of molecular gas at 10-103 K (e.g., Walter et al., 2002; Sturm et al., 2011; Bolatto et al., 2013). Most studies of galactic outflows focus on warm gas at $\approx$104K due to the strong emission and absorption lines in the rest-frame UV and optical. X-ray observations, however, can only be obtained for nearby objects, since only they are bright enough to be seen (e.g., Lehnert et al., 1999; Strickland & Heckman, 2009). This hot gas is fundamentally important as it dominates the thermal and kinetic energy of the outflow and provides the best evidence that winds actually escape their host galaxies (e.g., Lehnert et al., 1999).

Since galactic outflows are generally diffuse in nature and thus have low surface brightnesses in their extended emission, spatially-resolved emission line studies are difficult in anything other than nearby galaxies (e.g., Lehnert & Heckman, 1996; Westmoquette et al., 2009; Sharp & Bland-Hawthorn, 2010; Arribas et al., 2014; Spence et al., 2018). Therefore, absorption lines are the most suitable choice for studying the kinematics, column density, and mass and momentum outflow rates of winds.

Unfortunately, several issues complicate any analysis of the properties of winds. One of the persistent mysteries is the exact relationship between the hot and cold phases. One possibility is that cold material is being driven out of the host galaxy by ram pressure acceleration from hotter material entrainment (e.g., Lehnert & Heckman, 1996; Veilleux et al., 2005). However, this hypothesis runs into serious difficulties, both because: (i) shocks and conduction from the exterior medium tend to compress the cloud perpendicular to the direction of the flow, greatly reducing the momentum flux it receives; and (ii) instabilities and evaporation lead to rapid cloud disruption (Klein et al., 1994; Mac Low & Zahnle, 1994; Orlando et al., 2006, 2008; Scannapieco & Brüggen, 2015; Brüggen & Scannapieco, 2016). A second possibility is that cold gas is formed from the cooling of high temperature, $10^{7}-10^{8}$K gas, which is moving at high radial velocities. Such scenarios for explaining the cold gas in outflows was first proposed by Wang (1995a, b) and extended to include cooling (Efstathiou, 2000), infall and metallicity evolution (Silich et al., 2003, 2004; Tenorio-Tagle et al., 2007; Wünsch et al., 2011), and the formation of density inhomogeneities (Scannapieco, 2017). This possibility has received significant attention with recent observations of ultra-luminous infrared galaxies (ULIRGs) and other starbursts suggesting that it is plausible (Zhang et al., 2015; Martin et al., 2015; Thompson et al., 2016; Zhang et al., 2017).

Determining outflow properties without this full picture can lead to incorrect estimates in the mass outflow rate (Chisholm et al., 2018). For example, accurate rates can only be determined if we understand the ionization state of the outflowing gas. Photoionization models (e.g., Cloudy; Ferland et al., 2013) have been used interpret observations of UV absorption lines to estimate outflow rates. However, these models assume a steady state in which all the rates of ionizing and recombination are in equilibrium.

To better understand the impact of non-equilibrium effects, we present here a suite of one-dimensional spherically symmetric simulations that study the hydrodynamic and ionization structure of galactic outflows. We implement an inflow boundary condition in the simulation box that reproduces outflow properties of the freely expanding adiabatic outflow conditions at the sonic radius in the model of Chevalier & Clegg (1985). The ionization state of the gas is computed using a non-equilibrium atomic chemistry package that includes photoionization, photo-heating, and ion-by-ion recombination and cooling.

The structure of the paper is as follows: §2 introduces the galactic outflow model. §3 discusses the model framework and initial conditions. The results of our simulations are presented in §4. In §5, we discuss the N V-O VI ratio and its importance to observations. Finally, we summarize and conclude our study in §6.

2 Galactic Wind

We are interested in the ionization structure of an expanding outflow generated from a starbursting galaxy. Since the outflow moves over several to many gas scale heights per Myr and most starbursts last for many Myr (e.g., Greggio et al., 1998; Förster Schreiber et al., 2003; McQuinn et al., 2010), the outflow is often described as being in equilibrium. Such an assumption seems to match observations when reliable X-ray analyses can be made (e.g., Heckman et al., 1990, 1995; Ott et al., 2005; Strickland & Heckman, 2007; Yukita et al., 2012).

For a spherically symmetric outflow and assuming that gravitational forces are negligible the equations of motion become:

[TABLE]

[TABLE]

[TABLE]

where $r$ is the radial coordinate, $\rho$ is the mass density, $u$ is the radial velocity, $\gamma$ is the adiabatic index, and $P$ is the pressure.

The mass and energy input rates are:

[TABLE]

[TABLE]

where $R_{*}$ is the sonic radius (see below), $V=4\pi R_{*}^{3}/3$, $n_{i}$ is the number density of species $i$, $n_{e}$ is the electron number density, $\Lambda$ is total cooling function, and $\mathcal{H}$ is the total photo-heating rate.

In the case where $q_{e}$ and $q_{m}$ are zero outside the sonic radius the solution to these equations is presented in (Chevalier & Clegg, 1985, CC85 hereafter). These three parameters, the energy inflow rate, $\dot{E}$, the mass inflow rate, $\dot{M}$, and the sonic radius, $R_{*}$, define the hydrodynamics of the outflow. In the adiabatic case presented in CC85, the density, pressure, and velocity asymptotically approach $\rho\sim r^{-2}$, $P\sim r^{-10/3}$, and $u\sim u_{\infty}=\sqrt{2\dot{E}/\dot{M}}$. It is with this solution in mind that we perform the simulations presented here and implemented as the inflow boundary conditions in the simulation. In the following sections we describe our simulations in detail. In particular, the inflow boundary conditions meant to represent the outflow at the sonic radius.

3 Model Framework and Initial Conditions

All of the simulations presented here are run with MAIHEM (Gray et al., 2015; Gray & Scannapieco, 2016, 2017), our modified version of the adaptive mesh hydrodynamics code, FLASH(version 4.5 Fryxell et al., 2000). The atomic package has recently been updated to fully resolve some elements, i.e., carbon, nitrogen, oxygen, and neon, as well as now including argon. Our atomics package tracks the non-equilibrium evolution of 84 species across 13 elements: 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), iron (Fe I-Fe VI), and electrons (${e}^{-}$). Our atomic network consists of 310 reactions. For each species we consider electron impact ionization, radiative and dielectronic recombination, and photoionization due to a UV background.

To compute the photo-heating and photoionization rates, we utilize photoionization cross sections from Verner & Yakovlev (1995) for the inner shell electrons and Verner et al. (1996) for the outer shell electrons. The photoionization and photo-heating rates due to the UV background are parameterized by the line intensity at the Lyman limit, $\mathcal{J_{\nu}}=$10*-21J21* erg s*-1* cm*-2* Hz*-1* Sr*-1*, where J21 is a free parameter.

Two improvements have been made to our network since last presented in Gray & Scannapieco (2017). First, J21 is now defined cell-by-cell rather than as a global variable. This provides the ability to define a spatially varying, although temporally fixed, radiation field. Second, the photoionization and photoheating rates are computed at the beginning of every simulation. The radiation spectrum is defined at runtime via an input file. This allows for the flexibility of changing the input spectrum without recompiling the code.

The inclusion of the chemistry network and its constituent species necessitates updates to some of the equations presented above. The total cooling becomes

[TABLE]

where $\Lambda_{i}$ and $n_{i}$ are the cooling function and number density of species $i$, and $n_{e}$ is the number of electrons. The cooling rates are computed using Cloudy using the same procedure as described in Gray et al. (2015) and Gray & Scannapieco (2016), which has been updated to include the new ions in our expanded network. Our procedure creates a table of cooling rates that is then linearly interpolated over to compute the cooling rate for a given temperature. Similarly, the total photoionization rate, $\mathcal{H}$, is

[TABLE]

where $\mathcal{H}_{i}$ is the photoheating rate for ion $i$.

The mass inflow rate, $\dot{M}$, energy inflow rate, $\dot{E}$, and the sonic radius, $R_{*}$ define the properties of the galactic outflow and are used as the basis of our inflow boundary conditions. With these conditions, the energy flux of the outflow is then:

[TABLE]

where $v_{\rm outflow}$ is the outflow velocity and cs is the sound speed. At the sonic radius, R∗,

[TABLE]

The mass density of the outflow is,

[TABLE]

where $\Omega R_{*}^{2}$ is the effective area of the outflow. The density of the outflow is then given by, $\rho=\dot{M}/\Omega R_{*}^{2}v_{\rm outflow}$, for completely isotropic outflow, i.e., $\Omega=$4$\pi$. Here we take $\Omega=\pi$, as used in Scannapieco (2017) based on observations of a large sample of low-redshift starbursts (Heckman et al., 1990; Lehnert & Heckman, 1996; Martin et al., 2012).

The mass input rate can be scaled by the star-formation rate as, $\dot{M}=\beta\dot{M}_{SFR}$, where $\beta$ is the mass-loading factor of the outflow. Likewise, the energy input rate is scaled as $\dot{E}=\alpha\dot{E}_{SN}$, where $\alpha$ is the energy-loading factor, which accounts for the fraction of the energy from supernova directed into the outflow. If one assumes that a single supernova event generates 1051 ergs of energy, and that there is one supernova per 100 M⊙ formed, the total energy input rate is then $\dot{E}=$10${}^{49}\alpha\dot{M}_{SFR}$ erg yr*-1*.

The above outflow values are used to implement the boundary condition defined at the inner boundary. To complete the definition of all hydrodynamic variables, the temperature of the outflow is defined as,

[TABLE]

where mH is the mass of hydrogen, $\bar{A}$ is the average atomic weight, and kB is Boltzmann’s constant. The atomic composition of the gas inflowing at the boundary is assumed to have a solar composition. The initial ionization state of the gas is assumed to be equal to the collisional ionization equilibrium (CIE) values that depend only on the outflow temperature and the UV background radiation field. The inital CIE values depend on the initial hydrogen number density, input spectra, and ionization parameter. These values are model dependent and require the generation of a model specific CIE table. We employ Cloudy to compute these CIE values. This procedure creates a table of values that which is read in by FLASH at runtime and interpolated over to give the initial ionization state of the gas at the boundary.

3.1 Background UV Radiation Field

The incident UV background radiation field is a combination of two sources, the meta-galactic UV background from Haardt & Madau (2012) and a starburst model (SB99; Leitherer et al., 1999). The SB99 models are run with fixed stellar mass, with a metallicity of one-tenth solar, and the spectra that are evolved for $t=$4$\times$104 years. We have compared spectral energy distributions at earlier times and with solar metallicity and found little difference between them. Since the star formation rate is varied among our simulations, the SB99 spectral energy distributions is modified as,

[TABLE]

where $\dot{M}_{SFR}$ is the star-formation rate, tevo(=4$\times$104yr), is the time over which the spectrum is aged and the spectral energy distribution from SB99, $L_{6,{\rm SB99}}$, is normalized by a total stellar mass of 106 M*⊙. The conversion of SB99 spectral energy distributions with units of erg s-1* Å*-1* to spectral radiance is given as,

[TABLE]

where $\lambda$ is the wavelength, $c$ is the speed of light and R∗ is the sonic radius. Combining both the converted SB99 spectra and the meta-galactic background then creates the total composite spectrum. Figure 1 shows both components of the incident spectra as well as the set of composite spectra used in our simulations. As expected the star formation rate scales the spectra linearly in intensity. Overall, the meta-galactic background has a much lower intensity compared to the starburst background up to an energy of $\sim$100 eV. Although we use two components for the background spectra, we find that the photoionization and photoheating rates are dominated by the SB99 spectra and the meta-galactic background plays only a minor role.

The total background is spatially varied and follows an inverse square profile,

[TABLE]

where R*∗* is the sonic radius, Bν,SB99 is the Starburst99 component, and Bν,HM is the metagalactic background component. While the UV background does vary spatially, it is fixed in time. A minimum UV background is imposed that corresponds to a strength expected from the meta-galactic flux.

3.2 Simulation Setup and Nomenclature

Each simulation is run in one-dimensional spherical coordinates with the inner radius equal to the sonic radius, which is set to $R_{*}$=300 pc, consistent with M82 (e.g., McKeith et al., 1995; Strickland & Heckman, 2009). The maximum radius is 10 kpc. The ambient medium is set to an initial density that is much lower than the outflow density, $\rho_{\rm ambient}=\rho_{\rm outflow}/1000$ with an initial temperature of $T_{\rm ambient}=10^{6}$ K. The metallicity of the ambient medium is set to be equal to the outflow metallicity with ionization states set to their CIE values. These CIE values are computed assuming a negligible ionizing background. The actual values for the ambient medium are largely defined for completeness as the density of the ambient gas is low enough that it is removed by the outflow. The base grid of each simulation is comprised of 256 blocks. We allow up to four levels of refinement based on the density, pressure, and the radial velocity. This gives a maximum resolution of 2.38 pc. Each simulation is run until it reaches steady state, in which $v$, $\rho$, and $T$ remain constant.

Table 3 summarizes the simulations presented here. The mass inflow rate, the energy inflow rate, the sonic radius, metallicity, and whether or not the UV background is on are what define each simulation. The fiducial model is defined as a model with a mass inflow rate of $\dot{M}$=10 M*⊙* yr*-1*, an energy inflow rate of $\dot{E}$=1050 erg yr*-1*, and solar metallicity. For simplicity, each model is named by what is varied relative to our fiducial model. For example, Z05 represents a model with a mass inflow rate of $\dot{M}$=10 M*⊙* yr*-1*, an energy inflow rate of $\dot{E}$=1050 erg yr*-1*, and a gas metallicity of half that of the solar metallicity, and with the background UV field. The name is therefore Z05 because only the metallicity is changed relative to the fiducial model.

4 Results

We show the normalized density, pressure, and radial velocity for our fiducial model and the CC85 results (Figure 2). The CC85 results are obtained by using their equations (4) and (5) with the parameters of the fiducial model. The sound speed is given by, $c^{2}_{s}[M^{2}/2+1/(\gamma-1)]=q_{e}/q_{m}$, the outflow velocity is then $c_{s}M=u$, and the density is given by the constraint, $\rho ur^{2}=constant$. Following CC85, the density, pressure, and radial velocity are normalized by, $\rho^{*}=\rho/(\dot{M}^{3/2}\dot{E}^{-1/2}R_{*}^{-2})$, $P^{*}=P/(\dot{M}^{1/2}\dot{E}^{1/2}R_{*}^{-2})$, and $u^{*}=u/(\dot{M}^{-1/2}\dot{E}^{1/2})$ respectively.

We find very good agreement between the CC85 solution and our results. Note that since the inner radius of our simulations are placed at the sonic point, all of our results start at log10(R/R∗) =0. This comparison verifies that our inflow boundary conditions are accurate and that the density of our ambient medium is low enough that is does not substantially impede the outflow. For completeness, we show the gas temperature as a function of radius for each model (Figure 3). The initial temperature is given by equation (11).

4.1 Ionization state of the gas in the fiducial case

The utility of tracking the ionization state of individual ions is that we get an accurate picture of the evolution of the ionization state of the ensemble of elements in the outflow. Figure 4 shows the results for all ionization states of carbon. Here we compare the MAIHEM results with equilibrium results computed using Cloudy (Ferland et al., 2013), using our input ionizing background. Initially, Carbon is predominately fully ionized (i.e., C VII) and remains in a higher ionization state with respect to the equilibrium ionization at all radii (it is “over-ionized”). Near the gas injection boundary at the sonic radius, the densities are the highest, and reactions proceed quickly. Thus the relatively over-ionization is low. In this region, the equilibrium and full-chemistry results follow each other relatively closely. At larger radii, however, the recombination rates drop and non-equilibrium effects become more pronounced, often differing by orders of magnitude at distances above $\approx$10 times the sonic radius.

Note also that these effects are often difficult to discern given the partial information that is available from optical-UV spectroscopy of galaxies exhibiting outflows. For example at $R/R_{\star}\approx 10$, and C IV, C III, C II are all lower than their equilibrium values, and could be interpreted as corresponding to lower temperatures and/or lower ionization parameters without the additional information that the most of the carbon is in the form of C V and C IV rather than C I.

In Figure 5 we show the profile in atomic ionization states that are commonly observed. As the carbon figure, up to log10(R/R∗) $<$0.5, the outflow ionization states are well described by the equilibrium values. At larger radii, the full MAIHEM calculation generally produces higher ionization states than those predicted assuming collisional ionization equilibrium, but the fraction of a particular ion of an element is difficult to predict a priori. For example, N V, O VI, and Ne VIII all exceed their equilibrium ionization fractions and are over-ionized, while C IV, Si IV, and Mg II are much less over-ionized. These strong differences highlight the importance of considering non-equilibrium ionization when interpreting emission and absorption line spectra of galactic outflows and how outflows may cool and ultimately evolve.

4.2 The impact of varying the metallicity and UV background

We demonstrate the impact of varying the UV background on the ionization state of the outflowing gas in Figure 6. Only at large distances from the sonic radius, i.e.,  log10(R/R∗) $>$1, do some ionization states in the UV-free case deviate from the fiducial case. Higher ionization states, such as Ne VII, Si IV, and N V, show little difference between the fiducial and UV-free case and only for large radii. Lower ionization states, such as H I and Mg II, show large differences. These lower ionization states have larger photoionization cross sections, and therefore photoionization leads to higher ionization states. In the absence of a UV background field these low ionization states remain populated compared to the fiducial model.

A range of metallicities are found within galactic outflows. Numerical studies have found that these outflows can enrich the gas within galaxy halos up to supersolar metallicities (Choi et al., 2017) while others find the outflows are enriched up to 8Z*⊙* (Gan et al., 2018). Using a sample of quasars Xu et al. (2018) found that gas within the broad line region to be highly enriched. Conversely, outflows with low metallicities have been found in simulations of galaxy formation that employ realistic feedback perscriptions (e.g., Hafen et al., 2017; Muratov et al., 2017).

Our fiducial model assumes that the gas in the outflow has solar metallicity which we compare to a model with half solar metallicity and twice solar metallicity (Figure 7). Close to the sonic radius changing the metallicity has little effect on the resulting ionization state. Only at large distances from the sonic radius, log10(R/R∗) $>$0.75 do metallicity effects begin to impact the ionization state of the gas. In particular, lower metallicity leads to generally higher ionization states of certain elements compared to the fiducial model, *e.g.,*C III and C IV. Other elements, the reverse is true. Ionization states that are particularly impacted by lowering the metallicity having higher relative ionization are C III, C IV, N V, and Si IV.

In general, increasing the metallicity has the effect of creating a more neutral gas. For example, the Ne VIII and O VI curves show that these ions peak sooner and fall to much lower abundances at larger radii compared to the fiducial case. Others, such as Si IV and Mg II have similar profiles as the fiducial model, but have overall higher abundances. The impact of metallicity on altering the ionization state of the gas is due to the changes in recombination time. The recombination time scale is approximated as $\tau_{\rm rec}=1/Rn_{e}$ where $R$ is the total recombination rate which is a strong function of temperature and ne is the electron density. In general the recombination rate is smaller for hotter temperatures. Therefore, the recombination times are longer for Z05 because there are fewer electrons and the gas remains hotter. Z20, on the other hand, has shorter recombination times due to the increase in the electron densities and overall cooler gas temperatures, as shown in Figure 3.

4.3 The impact of varying the mass outflow rate

We show the effect of changing the mass inflow rate at the simulation boundary in Figure 8. Two models, M05 and M02, were calculated that decrease the mass inflow rate by 50 and 80% relative to the fiducial model while keeping the energy inflow rate constant (i.e., decreasing the relative amount of the mass-loading). Decreasing the mass injection rate has the effect of increasing the gas velocity that leads to a decrease in the gas density, increases the gas temperature, and increases the overall ionization state of the outflowing gas. For example, the relative fraction of the ion N V is reduced by two orders-of-magnitude in the M05 model compared to the fiducial model. This reduction in the fractional abundance of each state is seen across all species and is even more pronounced in the M02 model with its much lower mass injection rate. In fact, for M02 only the highest ionization states are found, such as N V, O VI  and Ne VII have any significant fractional abundance. The ionization states also have a much smoother profile in M05 and M02 compared to the fiducial model, that is, it does not have any peaks in ionization state found in the fiducial model. This is due to the higher initial outflow temperature which never falls below T$\approx$106K, which prevents the recombination to the ionization states we consider from becoming very efficient. Thus, the differences in ionization state are a simple consequence of longer recombination times for these atomic species. The outflow density in these models is lower than the fiducial model which decreases the number density of electrons, and consequently, increases the recombination time scale and the winds stay highly ionized.

4.4 The impact of varying the flow density

We next analyze the impact of varying both the mass and energy outflow rates while keeping the outflow velocity fixed (Table 3). This leads to outflows with varying initial mass densities but the same initial temperatures. Two models are run that vary both outflow rates (or inflow rates into the simulation). The model, M05E05, decreases both rates by a factor of two while the model, M50E50, increases them by a factor of five. This has the effect of raising or lowering the outflowing gas density while keeping the initial gas outflow velocity and temperature constant.

One additional change is made for M50E50 regarding the star formation rate. We note that these models are similar to those presented in Scannapieco (2017), where M05E05 represents conditions similar to the starburst galaxy, M82 (e.g., Förster Schreiber et al., 2003), while M50E50 represents conditions seen in ultra-luminous infrared galaxies (ULIRGs; Daddi et al., 2005; Piqueras López et al., 2012). In order to achieve such high mass outflow rates, the star-formation rate within the sonic radius must also be effectively increased. For M50E50, our input parameters are equivalent to assuming that the star-formation rate is a factor of ten higher than that of the M05E05 model. Due to the way we scale the UV background, increasing the star-formation rate also increasing the intensity of the ionizing radiation field leading to higher photoionization and photo-heating rates.

We show the results of these two models in Figure 9. At small radii, log10(R/R∗) $<$0.5, M05E05 matches very well with the fiducial model for the ionization states shown. At larger radii, M05E05 shows lower ionization mass fractions compared to the fiducial model. For example, both N V and C IV have peaks in their ionization that are slightly downwind of the fiducial model. Similar structure is seen in O VI and Ne VII but is more pronounced. In terms of the temperature, M05E05 is hotter at larger radii when compared to the fiducial model.

M50E50 shows much more dramatic differences with respect to the fiducial model than M05E05. In this model, the gas density is high enough that cooling and recombination become very efficient. In fact, at log10(R/R∗) $=0.5$ the gas temperature quickly decreases to T$\approx$104 K and leads to a unique distribution of ionization states where many of the higher ionization states quickly recombine (e.g., N V and Ne VIII) leading to high contributions to their elemental mass fractions for the ions, C III and Si IV. Interestingly, the contribution to the overall ionization of O from O VI is negligible for radii beyond log10(R/R∗) $\approx$0.5. At larger radii the ionization states remain largely unchanged, except for N V, suggesting that equilibrium values are reached. This equilibrium value is found in balancing the electron recombination rate to the photoionization rate. At T$\approx$104 K, the equilibrium temperature, electron impact ionization rates are typically small compared to either the electron recombination rates or the photoionization rates. As mentioned above, the recombination rate scales as the inverse of the electron number density and therefore has an $r^{-2}$ spatial profile, that is, the recombination timescale increases with increasing distance from the starburst. Similarly, the photoionization rate goes as $\Gamma_{i}=1/\gamma_{0,i}J_{21}$, where $\gamma_{0,i}$ is the normalized photoionization rate for ion $i$ and J21 quantifies the strength of the UV background. As such, it follows eq. (14) and has an $r^{-2}$ spatial profile. Therefore, both the recombination and photoionization time scales follow the same $r^{-2}$ profile and the ionization states are in equilibrium for log10(R/R∗) $>$0.5.

We also note that M50E50 may represent a case of catastrophic cooling within the outflow and may provide important insight into the development and evolution of superwinds and wind-driven superbubbles (e.g., Silich et al., 2007; Silich & Tenorio-Tagle, 2017). In particular, it may prove useful in understanding feedback processes in young super star clusters such as Mrk 71 (e.g., Oey et al., 2017). The interaction of the outflow with a hydrodynamically important ambient medium is outside the scope of this work and will form the basis of a future study.

4.5 Comparison to Schneider & Robertson (2018)

In the high-resolution galaxy outflow simulations presented in Schneider & Robertson (2018), the CC85 model is used as the basis for their outflow model.

Two outflow conditions are considered by these authors, a “high mass-loading” (SHigh) model with a mass and energy injection rate of $\dot{M}$=12 M*⊙* yr*-1* and $\dot{E}$=1.7$\times$1050 erg yr*-1*; and a “low mass-loading” (SLow) model with rates of $\dot{M}$=1.5 M*⊙* yr*-1* and $\dot{E}$=0.47$\times$1050 erg yr*-1*. As discussed by these authors, both cases represent physically motivated conditions. The “high-mass loading” model represents starburst-like conditions early in the burst where the star formation rate is high. The “low-mass loading” model represents a more mature starburst where the star formation rate is lower and gas in the ISM has been swept away. These two parameters set then span an interesting range in evolutionary time for a given starburst powered outflow.

We reproduced these two cases using our code (Figure 10). Compared to our fiducial model, the low mass-loading case produces higher ionization states at all radii. In fact, almost all of the studied ionization states are absent; only the highest ionization states are present but in very low abundance. This is reflected in the temperatures of each model that are much hotter in the simulations using the Figure 3 initial conditions than the fiducial model (Figure 3). These high temperatures create an environment where recombination rates are very low and the ionization state remains largely unchanged as the gas flows outward from the sonic radius.

The high mass-loading case, on the other hand, produces ionization states that are more similar to the fiducial case. This model produces peaks in N V, Ne VII, and O VI ionization states that are slightly higher radii from the sonic radius. For example, the peak of Ne VII is at log10(R/R∗) $\approx$0.75 in the fiducial model but found at log10(R/R∗) $\approx$1.0 in the high-mass loading case. This too is reflected in Figure 3 where the high-mass loading case remains slightly hotter than the fiducial model.

5 The column densities of N V and O VI

X-ray observations of galactic outflows show that the hot gas (T$>$107 K) dominates the kinetic and thermal energy of the outflow (e.g., Martin, 1999; Griffiths et al., 2000; Strickland & Stevens, 2000; Strickland & Heckman, 2009; Zhang et al., 2014). Since the hottest gas in outflows, $\sim$108 K, i.e., the “piston” that is driving the wind, is diffuse and has low emissivity, observations of this phase are particularly difficult. The best way to study the piston is to observe emission lines in the X-ray that are sensitive to gas at such high temperatures (Strickland & Heckman, 2009). Alternatively, one can study the strong UV/optical emission lines to probe the cooling rate and kinematics of outflowing gas at intermediate temperatures, 104 to 105.5 K (e.g., Lehnert & Heckman, 1996; Hayes et al., 2016). More commonly however, UV and optical absorption lines are used to probe the terminal velocities, column densities, covering fractions, and energetics of the outflowing gas over these intermediate temperatures (e.g., Heckman et al., 2015, 2017; Chisholm et al., 2018).

Chisholm et al. (2018) used the O VI 1032 Å doublet in the rest-frame far UV to study the column densities of warm/hot gas at temperatures $\approx$105.5 K and cooler. The authors used photoionization models to predict the abundance of O VI finding that photoionization alone was unable to explain the O VI column densities. But importantly for testing models of the gas ionization, they also found that N V is not detected in absorption even though N V 1243 Å is expected to be a relatively strong transition and probe similar gas. They found an upper limit on the N V-O VI column density ratio of $\approx$0.07.

As these photoionization models assume photoionization equilibrium, we are interested in the column densities of N V and O VI and their ratio produced by the suite of models presented here. The total column densities of our expanding wind models are mostly between about 7-24$\times$1014 and 7-33$\times$1013 cm*-2* for O VI and N V respectively. There are models with significantly lower column densities that are always those with very low mass outflow rates. Our estimated column densities are in good agreement with observations of nearby starburst and gas in the inner halos of galaxies (Grimes et al., 2009; Hayes et al., 2016; Chisholm et al., 2016; Bordoloi et al., 2016).

Figure 11 shows the ratio of the number densities of N V and O VI for each model. In order to prevent numerically unphysical ratios, only outflow regions where both species are found in appreciable amounts are considered as they are the regions that significantly contribute to the total column of these ions. A cutoff value of F${}_{i}>$10*-6* is used for the fractional ionization state of both ions (see the caption of Figure 5). At radii of log10(R/R∗) $<$0.75 nearly every model predicts the low abundance of N V compared to O VI. Beyond log10(R/R∗) $>$0.75, the presence of N V is model dependent. For example, the fiducial models shows the column density ratio of N V to O VI peaking above the upper limit of Chisholm et al. (2018) at log10(R/R∗) $>$0.9. Models with lower outflow densities and higher outflow velocities (and outflow temperatures), have ratios that consistent with this upper-limit, e.g., M05, M02, and SLow. The lack of a detection of a significant column of N V is a consequence of non-equilibrium effects and does not require alternative mechanisms of generating excess O VI. However, models with low densities generally do not agree with column densities observed in starbursts, as they are too low.

A plausible way that the models can be reconciled with observations, namely, producing sufficient O VI column densities while preserving a low column density ratio of N V and O VI is to increase the terminal velocity of the wind for constant mass outflow rate and/or have a moderately low outflow rate at constant wind velocity (models that appear roughly analogous to M82 - M05E05). In other words, within the context of our model, the mass loading of the wind is relatively modest. Even though some of the models are close to the upper limit observed in the column density (0.1 or less), a more complete set of models with more detailed physics and a wider range of observations of this ratio in starburst galaxies are necessary to see if this is quantitatively correct.

6 Summary and Conclusion

We have presented a suite of one-dimensional galaxy outflow models that track the non-equilibrium evolution of the ionization fractions of many species as a function of the sonic radius. The outflow model is based on the outflow model of Chevalier & Clegg (1985) and is implemented as a boundary condition. By using the Chevalier & Clegg (1985) model, the outflow is defined by three parameters at the boundary of the simulation, the mass inflow rate, $\dot{M}$, the energy inflow rate, $\dot{E}$, and the sonic radius, $R_{*}$. The inflowing gas is assumed to be in collisional ionization equilibrium as it enters the computational domain. In addition, we implement a spatially varying but constant photoionizing background based on Starburst99 models (Leitherer et al., 1999) and the meta-galactic UV background at z=0 (Haardt & Madau, 2012).

A range of models is presented that vary several parameters, including $\dot{M}$, $\dot{E}$, the metallicity of the gas, and the presence of the UV background. We find that every model reaches a steady state hydrodynamic solution that is well described by the normalized Chevalier & Clegg (1985) profiles of density, pressure, and outflow velocity. However, the distribution of ionization states of the species investigated is strongly dependent on the parameters of the outflow. Specifically, we find that reactions rates can not keep up with evolution of the outflowing material, and in many of the cases, the column densities of ions are much higher than one would predict assuming photoionization equilibrium.

To understand the non-equilibrium effects, we compared the ionization states from our fiducial model to collisional ionization equilibrium models from Cloudy (Ferland et al., 2013). Two regimes in the ionization are found. Near the outflows sonic radius, where the wind becomes freely expanding, the model ionization states are well described by equilibrium models. At a radius of log10(R/R∗) $\approx$0.75, which for the starburst galaxy, M82, is approximately 1.5 kpc (e.g., Strickland & Heckman, 2009, and references therein), the outflow gas has cooled enough that recombination becomes efficient for the higher ionization states. At larger radii, log10(R/R∗) $>$0.75, our models are generally have distributions in their ionization states that are skewed to higher ionization states compared to models that assume photoionization equilibrium. This is due to a combination of various non-equilibrium effects and but especially the long recombination times as the density decreases with radius.

One model in particular warrants further study. M50E50, characterized by large mass and energy inflow rates at the boundary, has the properties of an outflow undergoing catastrophic cooling. When compared to other models, M50E50 shows a drastic temperature falloff at log10(R/R∗) $\approx$0.6 (see Figure 3) and an ionization state distribution that seems to be in equilibrium with the ionizing background. Future work will study this phenomenon more closely and its relation to young super clusters (Oey et al., 2017) and interaction with a hydrodynamically important ambient medium.

Finally, we looked at the N V-O VI column density ratio. Chisholm et al. (2018) used the O VI 1032Å far UV absorption line to study the outflow mass loading in a $z\approx$ 2.9 galaxy. Interestingly, they were unable to detect the N V 1243Å absorption line even though it should be a prominent feature. The authors gave an upper limit of log(N V/O VI)$\approx$-1.2 and possible explanations for its absence. We have shown here that this ratio can be explained by the non-equilibrium nature of the ionization states within the outflow and is predicted over a wide range of outflow conditions.

The simulations presented here provide a key insight into the ionization state distribution within galactic outflows. Although the general structure of the ionization states is similar between models, the exact peaks in the abundance of a given ionization state of an element is model dependent. Thus, these models highlight the importance of accounting for non-equilibrium effects in studies of galaxy outflows.

We would like to thank Sally Oey and Todd Thompson for useful comments and discussions. Helpful comments by the referee are also gratefully acknowledged. The software used in this work was in part developed by the DOE NNSA-ASC OASCR Flash Center at the University of Chicago. E.S. was supported by NSF grants AST11-03608 and AST14-07835 and NASA theory grant NNX15AK82G. The figures and analysis presented here were created using the **yt ** analysis package (Turk et al., 2011).