Nonequilibrium Ionization States within Galactic Outflows: Explaining Their O vi and N v Column Densities

Galactic outflows are ubiquitous in intensely star-forming galaxies across all redshifts (Heckman et al. 1990, 2000, 2015; Lehnert & Heckman 1996; Pettini et al. 2002; Martin 2005; Rupke et al. 2005; Veilleux et al. 2005; Weiner et al. 2009; Bordoloi et al. 2014, 2016; Rubin et al. 2014; Chisholm et al. 2016). These outflows are powered by stars through supernovae (e.g., Mac Low & Ferrara 1999; Scannapieco et al. 2001, 2002; Mori 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; Agertz & Kravtsov 2015; Lu et al. 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 10⁷–10⁸ K material (Martin 1999; Strickland & Heckman 2007, 2009), near-UV and optical observations of warm ionized gas at ≈10⁴ (e.g., Pettini et al. 2001; Tremonti et al. 2007; Martin et al. 2012; Soto et al. 2012), and IR and submillimeter observations of molecular gas at 10–10³ 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 ≈10⁴ K 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⁷–10⁸ K gas, which is moving at high radial velocities. Such scenarios for explaining the cold gas in outflows was first proposed by Wang (1995a, 1995b) 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 (Martin et al. 2015; Zhang et al. 2017; Thompson et al. 2016).

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 to interpret observations of UV absorption lines to estimate outflow rates. However, these models assume a steady state in which all the rates of ionization and recombination are in equilibrium.

To better understand the impact of nonequilibrium 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 nonequilibrium atomic chemistry package that includes photoionization, photoheating, and ion-by-ion recombination and cooling.

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

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 megayear and most starbursts last for many megayears (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

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

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

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

where r is the radial coordinate, ρ is the mass density, u is the radial velocity, γ is the adiabatic index, and P is the pressure.

The mass and energy input rates are

$$\begin{eqnarray}{q}_{m}=\left\{\begin{array}{cc}\dot{M}/V & \mathrm{if}\ r\leqslant {R}_{* }\\ 0 & \mathrm{if}\ r\gt {R}_{* },\end{array}\right.\end{eqnarray} \tag{ 4 }$$

$$\begin{eqnarray}{q}_{e}=\left\{\begin{array}{cc}\dot{E}/V & \mathrm{if}\ r\leqslant {R}_{* }\\ -{\rm{\Lambda }}{n}^{2}-{ \mathcal H }n & \mathrm{if}\ r\gt {R}_{* },\end{array}\right.\end{eqnarray} \tag{ 5 }$$

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, Λ is total cooling function, and ${ \mathcal H }$ is the total photoheating rate.

In the case where q_e and q_m are zero outside of 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 ρ ∼ r⁻², P ∼ 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.

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 nonequilibrium 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 photoheating 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 photoheating rates due to the UV background are parameterized by the line intensity at the Lyman limit, ${{ \mathcal J }}_{\nu }$ = 10⁻²¹ J₂₁ erg s⁻¹ cm⁻² Hz⁻¹ Sr⁻¹, where J₂₁ is a free parameter.

Two improvements have been made to our network since last presented in Gray & Scannapieco (2017). First, J₂₁ 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

$$\begin{eqnarray}&&{\rm{\Lambda }}{n}^{2}=\displaystyle \sum _{i}{{\rm{\Lambda }}}_{i}{n}_{i}{n}_{e},\end{eqnarray} \tag{ 6 }$$

where Λ_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

$$\begin{eqnarray}&&{ \mathcal H }n=\displaystyle \sum _{i}{{ \mathcal H }}_{i}{n}_{i},\end{eqnarray} \tag{ 7 }$$

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:

$$\begin{eqnarray}&&\dot{E}=\dot{M}\left(1/2{v}_{\mathrm{outflow}}^{2}+3/2{c}_{s}^{2}\right),\end{eqnarray} \tag{ 8 }$$

where v_outflow is the outflow velocity and c_s is the sound speed. At the sonic radius, R_*,

$$\begin{eqnarray}&&{c}_{s}^{2}={v}_{\mathrm{outflow}}^{2}=\dot{E}/2\dot{M}.\end{eqnarray} \tag{ 9 }$$

The mass density of the outflow is

$$\begin{eqnarray}&&{\rm{\Omega }}{R}_{* }^{2}\rho {v}_{\mathrm{outflow}}=\dot{M},\end{eqnarray} \tag{ 10 }$$

where ${\rm{\Omega }}{R}_{* }^{2}$ is the effective area of the outflow. The density of the outflow is then given by $\rho =\dot{M}/{\rm{\Omega }}{R}_{* }^{2}{v}_{\mathrm{outflow}}$, for completely isotropic outflow, i.e., Ω = 4π. Here we take Ω = π, 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}}_{\mathrm{SFR}}$, where β is the mass-loading factor of the outflow. Likewise, the energy input rate is scaled as $\dot{E}=\alpha {\dot{E}}_{\mathrm{SN}}$, where α 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 10⁵¹ erg 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}}_{\mathrm{SFR}}$ erg yr⁻¹.

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

$$\begin{eqnarray}&&{T}_{\mathrm{outflow}}\approx {v}_{\mathrm{outflow}}^{2}{m}_{{\rm{H}}}\bar{A}/{k}_{{\rm{B}}},\end{eqnarray} \tag{ 11 }$$

where m_H is the mass of hydrogen, $\bar{A}$ is the average atomic weight, and k_B 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 initial 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 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 metagalactic 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 × 10⁶ yr. 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

$$\begin{eqnarray}&&{L}_{\mathrm{outflow}}={L}_{6,\mathrm{SB}99}\left(\displaystyle \frac{{\dot{M}}_{\mathrm{SFR}}\ast {t}_{\mathrm{evo}}}{{10}^{6}{M}_{\odot }}\right),\end{eqnarray} \tag{ 12 }$$

where ${\dot{M}}_{\mathrm{SFR}}$ is the star formation rate, t_evo(=4 × 10⁴ yr) is the time over which the spectrum is aged and the spectral energy distribution from SB99, and L_6,SB99 is normalized by a total stellar mass of 10⁶ M_⊙. The conversion of SB99 spectral energy distributions with units of erg s⁻¹ Å⁻¹ to spectral radiance is given as

$$\begin{eqnarray}&&{B}_{\nu ,\mathrm{SB}99}={L}_{\mathrm{outflow}}\left(\displaystyle \frac{{\lambda }^{2}}{c}\right)\left(\displaystyle \frac{1}{16{\pi }^{2}{R}_{* }^{2}}\right),\end{eqnarray} \tag{ 13 }$$

where λ is the wavelength, c is the speed of light, and R_* is the sonic radius. Combining both the converted SB99 spectra and the metagalactic 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 metagalactic background has a much lower intensity compared to the starburst background up to an energy of ∼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 metagalactic background plays only a minor role.

Zoom In Zoom Out Reset image size

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

$$\begin{eqnarray}&&{B}_{\nu }(R)=\left({B}_{\nu ,\mathrm{SB}99}+{B}_{\nu ,\mathrm{HM}}\right){\left(\displaystyle \frac{{R}_{* }}{R}\right)}^{2},\end{eqnarray} \tag{ 14 }$$

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 metagalactic 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, ρ_ambient = ρ_outflow/1000 with an initial temperature of T_ambient = 10⁶ 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, ρ, and T remain constant.

Table 1 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⁻¹, an energy inflow rate of $\dot{E}$ = 10⁵⁰ erg yr⁻¹, 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⁻¹, an energy inflow rate of $\dot{E}$ = 10⁵⁰ erg yr⁻¹, 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.

Table 1.  Summary of the Simulations Presented

Name	$\dot{M}$	$\dot{E}$	α	β	Voutflow	n	T	Z	UV	SFR
	(M⊙ yr−1)	(1050 erg yr−1)			(km s−1)	(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

Note. The first column gives the name for each model. The second and third columns give the mass inflow rate and energy inflow rate respectively. Columns 4 and 5 give the energy loading and mass loading constants. Columns 6–8 present the resulting outflow velocity, number density, and temperature of the outflow. The ninth column gives the metallicity of the outflow. A checkmark in column 10 shows a model with the UV background implemented. Finally, column 11 gives the star formation rate, which is used in determining the magnitude of the UV background, see Equation (12).

Download table as:  ASCIITypeset image

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}_{s}^{2}[{M}^{2}/2+1/(\gamma -1)]={q}_{e}/{q}_{m}$, the outflow velocity is then c_sM = u, and the density is given by the constraint, ρur² = 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.

Zoom In Zoom Out Reset image size

We find very good agreement between the CC85 solution and our results. Note that since the inner radius of our simulations is placed at the sonic point, all of our results start at ${\mathrm{log}}_{10}(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 it 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).

Zoom In Zoom Out Reset image size

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 "overionized"). Near the gas injection boundary at the sonic radius, the densities are the highest, and reactions proceed quickly. Thus the relatively overionization 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 nonequilibrium effects become more pronounced, often differing by orders of magnitude at distances above ≈10 times the sonic radius.

Zoom In Zoom Out Reset image size

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_⋆ ≈ 10, and C iv, C iii, and 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 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 ${\mathrm{log}}_{10}(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 CIE, 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 overionized, while C iv, Si iv, and Mg ii are much less overionized. These strong differences highlight the importance of considering nonequilibrium ionization when interpreting emission and absorption line spectra of galactic outflows and how outflows may cool and ultimately evolve.

Zoom In Zoom Out Reset image size

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., ${\mathrm{log}}_{10}(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.

Zoom In Zoom Out Reset image size

A range of metallicities is 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 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 prescriptions (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, ${\mathrm{log}}_{10}(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. For 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.

Zoom In Zoom Out Reset image size

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 timescale is approximated as τ_rec = 1/Rn_e where R is the total recombination rate which is a strong function of temperature and n_e 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, an increase in the gas temperature, and an increase in 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 that never falls below T ≈ 10⁶ K, 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 timescale, and the winds stay highly ionized.

Zoom In Zoom Out Reset image size

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 1). 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 in 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 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 10 higher than that of the M05E05 model. Due to the way we scale the UV background, increasing the star formation rate also increases the intensity of the ionizing radiation field leading to higher photoionization and photoheating rates.

We show the results of these two models in Figure 9. At small radii, ${\mathrm{log}}_{10}(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.

Zoom In Zoom Out Reset image size

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 ${\mathrm{log}}_{10}(R/{R}_{* })$ = 0.5 the gas temperature quickly decreases to T ≈ 10⁴ 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 ${\mathrm{log}}_{10}(R/{R}_{* })$ ≈ 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 ≈ 10⁴ 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⁻² spatial profile; that is, the recombination timescale increases with increasing distance from the starburst. Similarly, the photoionization rate goes as Γ_i = 1/γ_0,iJ₂₁, where γ_0,i is the normalized photoionization rate for ion i and J₂₁ quantifies the strength of the UV background. As such, it follows Equation (14) and has an r⁻² spatial profile. Therefore, both the recombination and photoionization timescales follow the same r⁻² profile and the ionization states are in equilibrium for ${\mathrm{log}}_{10}(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⁻¹ and $\dot{E}$ = 1.7 × 10⁵⁰ erg yr⁻¹; and a "low mass-loading" (SLow) model with rates of $\dot{M}$ = 1.5 M_⊙ yr⁻¹ and $\dot{E}$ = 0.47 × 10⁵⁰ erg yr⁻¹. 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.

Zoom In Zoom Out Reset image size

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 ${\mathrm{log}}_{10}(R/{R}_{* })$ ≈ 0.75 in the fiducial model but found at ${\mathrm{log}}_{10}(R/{R}_{* })$ ≈ 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.

X-ray observations of galactic outflows show that the hot gas (T > 10⁷ 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, ∼10⁸ 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, 10⁴–10^5.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 ≈10^5.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 ≈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) × 10¹⁴ and (7–33) × 10¹³ cm⁻² 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; Bordoloi et al. 2016; Chisholm et al. 2016; Hayes 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⁻⁶ is used for the fractional ionization state of both ions (see the caption of Figure 5). At radii of ${\mathrm{log}}_{10}(R/{R}_{* })$ < 0.75 nearly every model predicts the low abundance of N v compared to O vi. Beyond ${\mathrm{log}}_{10}(R/{R}_{* })$ > 0.75, the presence of N v is model dependent. For example, the fiducial model shows the column density ratio of N v to O vi peaking above the upper limit of Chisholm et al. (2018) at ${\mathrm{log}}_{10}(R/{R}_{* })$ > 0.9. Models with lower outflow densities and higher outflow velocities (and outflow temperatures), have ratios that are 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 nonequilibrium 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.

Zoom In Zoom Out Reset image size

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.

We have presented a suite of one-dimensional galaxy outflow models that track the nonequilibrium 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 CIE 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 metagalactic 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 cannot 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 nonequilibrium effects, we compared the ionization states from our fiducial model to CIE models from Cloudy (Ferland et al. 2013). Two regimes in the ionization are found. Near the outflow's sonic radius, where the wind becomes freely expanding, the model ionization states are well described by equilibrium models. At a radius of ${\mathrm{log}}_{10}(R/{R}_{* })$ ≈ 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, ${\mathrm{log}}_{10}(R/{R}_{* })$ > 0.75, our models 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 nonequilibrium effects and 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 ${\mathrm{log}}_{10}(R/{R}_{* })$ ≈ 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 ≈ 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) ≈ −1.2 and possible explanations for its absence. We have shown here that this ratio can be explained by the nonequilibrium 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 nonequilibrium 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).

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