Simulating Gas Inflow at the Disk-Halo Interface

TL;DR

This study uses high-resolution hydrodynamic simulations to explore how inflowing gas interacts with galactic outflows at the disk-halo interface, revealing cold cloud formation, partial accretion, and observable signatures consistent with the Milky Way.

Contribution

It provides new insights into the fragmentation, dynamics, and observational signatures of gas inflow interacting with outflows at the disk-halo boundary using detailed simulations.

Findings

Cold gas clouds form from inflow due to cooling and instabilities.

Injected gas can partially accrete onto the disk, with efficiency depending on initial velocity.

Simulated absorption/emission signatures match observed properties of the Milky Way.

Abstract

The interaction between inflowing gas clouds and galactic outflows at the interface where the galactic disk transitions into the circumgalactic medium is an important process in galaxy fueling, yet remains poorly understood. Using a series of tall-box hydrodynamic ENZO simulations, we have studied the interaction between smooth gas inflow and supernovae-driven outflow at the disk-halo interface with pc-scale resolution. A realistic wind of outflowing material is generated by supernovae explosions in the disk, while inflowing gas is injected at the top boundary of the simulation box with an injection velocity ranging from $10-100 \rm \ km \ s^{-1}$. We find that cooling and hydrodynamic instabilities drive the injected gas to fragment into cold ($\sim 10^{3}$ K) cloud clumps with typical densities of $\sim 1 \rm \ cm^{-3}$. These clumps initially accelerate before interacting and partially mixing with the outflow and decelerating to velocities in the 50-100 $\rm km \ s^{-1}$ range. When the gas clumps hit the disk, $10\%-50 \%$ of the injected material is able to accrete (depending on the injection velocity). Clumps originating from gas injected with a higher initial velocity approach the disk with greater ram pressure, allowing them to penetrate through the disk in low density regions. We use (equilibrium) CLOUDY photoionization models to generate absorption and emission signatures of gas accretion, finding that our mock HI and H$\alpha$ observables are prominent and generally consistent with measurements in the Milky Way. We do not predict enhanced emission/absorption for higher ionization states such as OVI.

Figures

Tables

Table 1. Table 1: Summary of gas injection parameters in each simulation.

Name	Velocity [km s-1]	Density [${\mathrm{cm}}^{-3}$]	Temperature [K]	Start Time [Myr]	Stop Time [Myr]	Mass† [g]
BURST50	50	0.05	${10}^4$	50	75	$2.10\times {10}^5{M}_{\odot }$
BURST100	100	0.05	${10}^4$	50	62	$2.11\times {10}^5{M}_{\odot }$
IMD10	10	0.05	${10}^4$	0	70	$1.04\times {10}^5{M}_{\odot }$

Equations

$$n_{Xi}=n_H(n_X/n_H)(n_{Xi}/n_X)[\mathrm{c}{\mathrm{m}}^{-3}]$$

$$\Gamma =\frac{ϵ(\rho ,T)}{4\pi }\frac{\lambda }{hc}[\mathrm{photons}{\mathrm{s}}^{-1}\mathrm{c}{\mathrm{m}}^{-2}\mathrm{s}{\mathrm{r}}^{-1}]$$

Full text

galaxies: evolution — ISM — hydrodynamics — methods: numerical

1 Introduction

Gas accretion is needed throughout the cosmic evolution of spiral galaxies in order to replenish the galactic disk with hydrogen and fuel long-term star formation (Erb, 2008; Hopkins et al., 2008; Putman et al., 2009). Without a fresh influx of low metallicity gas (Chiappini et al., 2001; Fenner & Gibson, 2003; Sommer-Larsen et al., 2003), most star-forming galaxies including the Milky Way will exhaust their supply of disk gas (Chomiuk & Povich, 2011; Leroy et al., 2013), although stellar outflows can delay this (Leitner & Kravtsov, 2011).

There are a variety of external gas reservoirs beyond the galactic disk that may contain ample gas for sustaining a constant star formation rate. The majority of low-redshift baryons lie in the intergalactic filaments that make up the large-scale structure (Cen & Ostriker, 1999; Shull et al., 2012). These filaments feed into the potential wells around galaxies, where circumgalactic gas pools to form the galactic halo. Due to its proximity to the galactic disk, this circumgalactic medium (CGM) is a direct source of fuel for future star formation and closely influences galaxy evolution (Putman et al., 2012; Tumlinson et al., 2017). The CGM is a multi-phase combination of inflowing material from the intergalactic medium (IGM), gas stripped from orbiting satellite clouds, and enriched outflow expelled from the galactic disk. All of these sources are thought to fuel the galactic disk but with differing accretion processes (Putman, 2017).

Cosmological simulations have demonstrated that gas enters into the CGM along the cosmic web, resulting in filamentary inflow within the hotter, more diffuse halo gas. This inflow can be primarily cold when the halo virial temperature is relatively low (and so cooling is effective) and/or at high redshift, when the filament density is higher (e.g., Kereš et al., 2005). The inflow transitions to hotter gas flows at higher masses and lower redshifts. There is evidence from recent simulations that such flows heat, either via shocks or mixing before entering the galaxy (Nelson et al., 2013), although the numerical resolution in such cosmological simulations tends to be relatively poor in the CGM. Cosmological numerical models of halos with masses similar to the Milky Way suggest that low-redshift accretion is dominated by warm-hot ($10^{5}<$ T $<10^{6}$ K) filamentary inflow at $50-150$ km s*-1* (Joung et al., 2012b); such inflow fragments and mixes with feedback from the galactic disk (Joung et al., 2012b; Fernández et al., 2012).

The majority of inflowing gas is expected to accrete towards the outskirts of a galaxy, avoiding strong feedback form the galactic center (Stewart et al., 2011; Fernández et al., 2012). Even towards the outskirts of a galaxy, gas will encounter stellar feedback on its path to the disk. In cosmological hydrodynamic simulations, interactions between filamentary inflow and the surrounding halo/outflowing gas are important for creating over-densities in the filaments that facilitate gas cooling and fragmentation (Kereš et al., 2009; Joung et al., 2012b; Fernández et al., 2012). Interactions between inflowing and outflowing gas at the disk-halo interface are likely similarly important in the final stages of gas accretion, but are difficult to study with the kpc-scale resolution in cosmological zoom-in hydrodynamic simulations.

All gas exchanged between the galaxy and the surrounding halo must pass through a boundary layer where the galactic disk transitions into the lower galactic halo. This region is called the disk-halo interface (Norman & Ikeuchi, 1989; Ford et al., 2010; Putman et al., 2012), and it is a crucial boundary for studying the final stages of gas accretion. Across the disk-halo interface, there is a transition in the phase, density, and kinematic properties of the gas. A layer of dense, cold HI gas makes up most of the disk mass in the solar neighborhood, ranging in thickness from $<100$ pc in the galactic center to $\sim 1$ kpc along the Galactic Warp (Dickey & Lockman, 1990; Levine et al., 2006). This neutral component transitions into a diffuse warm ionized component (or Reynolds Layer) extending to a height of $\sim 2$ kpc (Reynolds, 1993; Haffner et al., 2003; Gaensler et al., 2008). An extensive population of cold HI clouds have been found embedded in this warm gas layer within $1-2$ kpc of the galactic disk (Lockman, 2002; Ford et al., 2010; Peek et al., 2011; Saul et al., 2012). The volume above the diffuse warm ionized layer is filled by warm-hot/hot gas blending into the Galactic halo (Miller & Bregman, 2013, 2015)

The gas distribution at the disk-halo interface is likely to be quite dynamic, due to multiple effects. This includes the impact of feedback expelling hot gas to large scale heights as well as possible gas cooling close to the galactic disk (Joung & Mac Low, 2006; Hill et al., 2012; Putman et al., 2012; Creasey et al., 2013; Li et al., 2017a). Additional evidence of gas cooling is provided by the observed HI cloud populations associated with the disk-halo interface (Lockman, 2002; Ford et al., 2010; Peek et al., 2011; Saul et al., 2012). These clouds are thought to condense out of the surrounding medium, and a variety of models have been proposed to facilitate this gas cooling (Bregman, 1980; Norman & Ikeuchi, 1989; Houck & Bregman, 1990; Fraternali & Binney, 2008; Marasco et al., 2012; Joung et al., 2012a). The survival of these clouds, and therefore their ability to fuel the galactic disk, is dependent on their mass, their height above the disk, and their interactions with the surrounding gas as well as galactic outflow (Heitsch & Putman, 2009).

Extragalactic observations around spiral galaxies reveal multiple pieces of evidence supporting the presence of gas inflow. One of the strongest pieces of direct evidence for accretion comes from absorption line spectra of UV-bright stars across the disk of M33. Since M33 is inclined at $55^{\circ}$, these observations are able directly probe the vertical gas kinematics, revealing disk-wide ionized gas inflow (Zheng et al., 2017). At greater distances, the disk-halo interface has been observed in handful of edge-on extragalactic systems. These observations probe the vertical velocity structure of the gas indirectly through measurements of the radial velocity profile. Gas in the halos of edge-on galaxies exhibits a kinematic transition where both neutral gas (Sancisi et al., 2001; Fraternali et al., 2002; Oosterloo et al., 2007; Heald et al., 2011) and ionized gas (Rand, 2000; Heald et al., 2006, 2007; Bizyaev et al., 2017) decrease in rotational velocity as a function of distance from the galactic disk. This rotational velocity gradient suggests the presence of gas inflow at the disk-halo interface, as inflowing gas should lag the rotational velocity of the disk and accreting gas merging with the galactic disk should be co-rotating. In the Milky Way, this trend is observed in the the populations of discrete intermediate velocity HI clouds observed at the disk-halo interface, which are known to exhibit net inflow velocities (Wakker, 2001; Lockman, 2002; Ford et al., 2010; Peek et al., 2011; Saul et al., 2012). In a recent sample, these HI clouds can be split into a kinematically warm population lagging the rotation of the disk at larger scale heights and a kinematically cold co-rotating population within 1 kpc (Saul et al., 2012).

The presence of outflowing gas at the disk-halo interface regulates gas accretion and possibly even re-fuels the galactic disk in a fountain cycle of expulsion and re-accretion (Shapiro & Field, 1976; Bregman, 1980). Multiple feedback mechanisms expel gas from the galactic disk, influencing the evolution of the galaxy as well as the composition of the circumgalactic medium. Radiative feedback sends ionizing photons into the halo while mechanical feedback adds material, energy, and momentum from the disk back into the surrounding CGM. One of the most energetic and ubiquitous forms of mechanical feedback is stellar feedback in the form of overlapping supernovae explosions (Joung & Mac Low, 2006; Creasey et al., 2013; Fielding et al., 2018), which can launch multi-phase outflows of metal-rich gas into the CGM with velocities of $\sim 1-100$ km s*-1* (Li et al., 2017a; Kim et al., 2017).

We present a suite of small box hydrodynamic ENZO simulations designed to study the interaction between inflowing and outflowing gas at the disk-halo interface with pc-scale resolution. This work builds upon simulations of gas outflow presented in Li et al. (2017a), where they studied the energy, mass, and metallicity of gas blown out of the galactic disk by overlapping supernovae explosions. We have added gas inflow to these simulations, the setup for which is described in Section 2. Section 3 presents the results of our simulation suite, and an overview of the gas evolution, cloud fragmentation, and eventual accretion outcomes. We use this simulation suite in conjunction with photoionization modeling to make observational predictions in Section 4. Finally, we discuss the implications of this work and and summarize our findings in Sections 5 and 6 respectively.

2 Methods

The simulations in this paper were performed using the adaptive-mesh hydrodynamics code Enzo (Bryan et al., 2014). Each simulation is a small box centered on the galactic plane with supernovae-driven outflow throughout the simulation volume. The initial setup as well as the supernova outflow prescription are identical to the $\rm\Sigma 10-KS$ model detailed in Li et al. (2017a). We refer to that paper for details, but also briefly describe the initial setup in Section 2.1 below. In this work we have augmented that work by injecting inflowing gas at the top boundary of the box in order to study the interaction between inflowing and outflowing gas at the disk-halo interface; inflow characteristics are described in more detail in Section 2.2.

2.1 Initial Setup

These simulations focus on a local volume of gas centered on the galactic disk. Each simulation is initialized as a rectangular box (elongated along the z-axis) with [x,y,z] dimensions $[0.35\times 0.35\times 5]$ kpc. A disk of gas is set up in the x-y plane at the center of the box (z = 0), and a stratified density halo extends $\pm 2.5$ kpc above and below the disk (with an initial temperature of $10^{4}$ K). The x and y boundaries of the box remain periodic throughout the simulation, and the z-axis boundaries are initially set to outflow. After a time delay (given in Table 1 under “Start Time”), the top z-axis boundary is changed to inflow and we inject a 1 kpc-long filament of gas inflowing towards the galactic disk. We use adaptive mesh refinement (AMR) to resolve the detailed interactions between the inflowing and outflowing gas. There are two levels of grid refinement within the simulation volume, one at $\pm 1$ kpc from the disk and another at $\pm 0.5$ kpc. Each sub-grid increases the resolution by a factor of two. The finest grid within $\pm 0.5$ kpc of the galactic disk achieves a spatial resolution of 2 pc.

We model the galactic disk near the solar neighborhood with a gas surface density ($\rm\Sigma_{gas}$) of $\rm 10\ M_{\odot}\ pc^{-2}$, corresponding to a SFR surface density ($\rm\dot{\Sigma}_{SFR}$) of $6.31\times 10^{-3}\rm\ {M}_{\odot}^{-1}$  kpc*-2* yr*-1* on the Kennicutt-Schmidt (KS) relation (Kennicutt, 1989). The photoelectric heating rate is set to a constant value of $\rm 1.4\times 10^{-26}\ ergs\ s^{-1}$ per H atom. The gas is initially isothermal and distributed in hydrostatic equilibrium, with gravity modeled as the combination of a baryonic stellar disk and an NFW dark matter halo. We apply a gas cooling curve appropriate for half solar-metallicity gas, ranging in temperature from $300-10^{9}$ K (Rosen & Bregman, 1995; Tasker & Bryan, 2006).

In order to drive inhomogeneous outflows, supernovae are injected near the midplane. Each supernova adds $10\ M_{\odot}$ of gas and $10^{51}$ ergs of thermal energy. The horizontal (x-y) distribution of supernovae is random. We implement both core-collapse and Type Ia supernovae using two different vertical distributions. Most of the supernovae ($90\%)$ are core-collapse, modeled as a gaussian vertical (z) distribution with a scale height $\rm h_{cc}=150$ pc. The remaining $10\%$ are Type Ia supernovae distributed exponentially along the z-axis with a scale height $\rm h_{Ia}=325$ pc. Supernovae are set off at a constant rate of $\dot{\Sigma}_{SFR}/m_{0}=\rm 4.2\times 10^{-5}\ kpc^{-2}\ yr^{-1}$, where we assume one SN explosion for every $m_{0}=150\ M_{\odot}$ of star formation. This supernovae prescription drives material towards the outflow z-axis boundaries of the simulation box. See Li et al. (2017a) for more details of the initial setup as well as a detailed analysis of the resulting outflows.

2.2 Inflow Conditions

In order to inject a cloud into the domain, the top z-axis boundary condition is changed from outflow to inflow. This is generally done after a time delay (see Table 1) in order to allow the outflow to become established. Supernovae continue to expel gas towards the boundaries, but this outflow is now met by a layer of inflowing material. All injected gas at the boundary is given an initial density, temperature, and z-velocity towards the disk. The initial values for these injected gas parameters are set by a number of considerations, including the desired inflow mass flux as well as the ram pressure of the outflowing gas. We characterize the inflow rate in terms of the star formation rate surface density, motivated by the idea that the inflowing gas should more than replenish the ongoing star formation in the disk. Accordingly, we compare the injected mass flux ($\rho v$) with the star formation rate surface density ($\rm\dot{\Sigma}_{SFR}$). For example, if inflows of a given rate had a (temporally and spatially averaged) covering fraction of $10\%$, then the inflow rate (when “on”) required to maintain star formation, on average, would be $10$ times the mean star formation rate surface density. Based on this logic, we explore ranges from $2-20$ times the mean star formation rate surface density. An additional constraint is that the ram pressure of the injected gas must be greater than the ram pressure of the outflowing gas. If this condition is not met, the injected gas will be unable to enter the box, an issue which we discuss more in the following explanation of Figure 1.

Figure 1 (top) shows the combinations of gas density and z-velocity resulting in an injected mass flux equivalent to $2\times\rm\dot{\Sigma}_{SFR}$ (solid), $10\times\rm\dot{\Sigma}_{SFR}$ (dashed), and $20\times\rm\dot{\Sigma}_{SFR}$ (dotted). The orange circle, purple square, and pink triangle denote the initial values used in each of our three simulation runs (See Table1); they fall on the three given mass flux isopleths respectively. Figure 1 (bottom) shows 1D profiles of the ram pressure (black solid) and thermal pressure (black dashed) of the outflowing gas at 50 Myr. By this time, the supernovae explosions have set up a continuous (steady-state) stream of inhomogeneous outflow at the z-axis boundaries. At these boundaries, the ram pressure ($\rm P_{ram}=\rho v^{2}$) of the steady-state outflowing gas dominates the total pressure ($\rm P_{ram}+P_{thermal}$). The ram pressure of the injected gas for each simulation is shown in solid colors corresponding to the legend in Figure 1a. The ram pressure of the 50 km s*-1* and 100 km s*-1* injected gas exceeds the ram pressure of the outflowing gas, allowing injected gas with these initial conditions to enter the box at 50 Myr. Gas injected at 10 km s*-1* comes close to equating the mass flux with the star formation rate surface density ($\rho v=\rm\dot{\Sigma}_{SFR}$), but the ram pressure of this low-velocity injected gas falls a factor of 3 below the ram pressure of the outflowing gas. As described in Section 2.3, we manage to simulate this low-velocity inflow by immediately injecting gas into the simulation box at the first time step, when the outflow ram pressure is much lower.

In addition to setting the initial density, temperature, and z-velocity of the injected gas, we have also added a colour density field to tag the injected gas and track its progression throughout the simulation volume. The colour density is the density of only the injected gas (from the top boundary), while the rest of the gas in the domain starts with a near-zero colour density field (including SN injected gas in the mid-plane). We can divide the colour density by the total gas density to create a normalized colour density field, which represents the fraction of inflowing gas present in a cell volume. In our analysis we use both the colour density and the normalized colour density. In some analysis, we use a normalized colour density cut of $>50\%$ to separate the properties of the injected gas from those of the outflowing gas.

2.3 Simulation Overview

We have run three variations of our basic simulation: BURST50, BURST100, and IMD10. Each follows the initial setup described in Section 2.1 but with different inflow rates. The number in the naming scheme for each simulation is the z-velocity of the injected gas in km s*-1*. The injected gas is given an initial density (n = 0.5 $\rm cm^{-3}$) and temperature (T = $10^{4}$ K), both of which are kept constant across all runs. The runs are differentiated by the start time of the gas injection, the z-velocity of the injected gas, and the injection stop time. The injection duration (stop $-$ start) is set so that the inflowing gas in each simulation reaches a length of $\sim 1$ kpc, corresponding approximately to a mass of $10^{5}$ M⊙. We halt the gas injection by reverting the top z-boundary condition to outflow, creating a finite burst of injected gas inflowing towards the galactic disk.

The BURST50 and BURST100 simulations inject a burst of inflowing gas at 50 km s*-1* and 100 km s*-1* respectively. Both of these runs simulate inhomogeneous supernovae-driven outflows interacting with an inflowing $\sim 1$ kpc-long burst of injected gas. We delay the injection start time by 50 Myr in order to set off enough supernovae to develop a steady-state outflow at the box boundaries. The IMD10 run immediately injects gas at 10 km s*-1* starting from the first time step. We are able to lower the inflow velocity in this run because steady state outflow has not yet been established at early times, the outflow ram pressure is much lower, and the 10 km s*-1* injected gas enters the box relatively unopposed. The injected gas is initially only in contact with the halo and later impacts the outflowing gas closer to the disk. We again stop injecting gas when the 10 km s*-1* inflow reaches about 1 kpc in length.

3 Results

We present the results of our three primary simulations: the BURST50, BURST100, and IMD10 runs. This analysis was performed using the python package yt (Turk et al., 2011). We begin with an overview of the physical gas structure (Section 3.1) and the time evolution throughout each simulation (Section 3.2). We then detail the fragmentation of the injected gas (Section 3.3) and its ultimate fate upon reaching the disk (Section 3.4).

3.1 Physical Structure

We begin by briefly describing how the inflowing, high density gas reacts as it impacts the outflowing, lower-density material. This situation produces the classic three-wave structure familiar from supernovae studies. When we inject gas at the top of the simulation box, we create a contact discontinuity at the boundary between the outflowing and injected material, as well as a weak “forward” shock propagating down into the outflowing gas, and a relatively high-density “reverse” shock of decelerated injected gas, which piles up behind the cloud front. This reverse shock propagates up into the inflowing wind (which is assumed to be relatively cool) but at a speed lower than the inflow rate so that the reverse shock is swept down towards the galaxy.

This structure can initially be clearly seen at the upper boundary condition immediately after injection (not shown), but rapidly evolves into a more complicated mixture of clouds as instabilities set in. As is well known, higher density gas above lower density gas at a contact discontinuity in a gravitational field creates a boundary that is Rayleigh-Taylor (RT) unstable. The reverse shock in the injected gas amplifies the growth of the RT instability by increasing the density disparity at the contact discontinuity. Rayleigh-Taylor fingers form from the cold, high-density injected gas. Over time, these features grow in length, and shearing with the outflowing gas triggers the Kelvin-Helmholtz (KH) instability. The RT and KH instabilities in combination lead to the eventual fragmentation of the injected gas into smaller clumps (Section 3.3) that will either mix into the outflowing gas, pass through the disk, or successfully accrete (Section 3.4). This rapid cloud formation means that the precise form of injected gas is unlikely to dramatically change the results we find below.

3.2 Global Evolution

The details of the global gas evolution are explored in Figures 2$-$4. Figure 2 shows slices in density, temperature, pressure, z-velocity, colour density and normalized colour density for four different time steps throughout the BURST50 simulation. At the boundary where the injected gas impacts the outflow, the injected gas decelerates rapidly and the z-velocity slices show a layer of stalled material. This corresponds to an over-dense layer in the injected gas density exceeding the initial conditions. In the temperature slices, this density enhancement allows the gas layer to cool efficiently, driving the temperature to $10^{4}$K and below. Rayleigh-Taylor fingers of cold, over-dense injected gas are seen falling into the less-dense halo environment at all stages of the time evolution in Figure 2. These inflowing fingers fill low-density voids in the inhomogeneous outflow where there is only halo gas and the resulting shear drives Kelvin-Helmholtz instabilities with the adjacent outflowing gas, resulting in smaller fragments and generating strands of small clouds. The colour density and normalized colour density (lower two sets of panels) track the progress of the inflowing cloud as it fragments and mixes with the outflowing gas. We discuss the properties of these fragmented cloud structures in more detail in Section 3.3.

The density and temperature slices in Figure 2 show the multiphase structure generated by this interaction. The corresponding pressure slices generally show less variation on small scales (although by no means none: a rich set of shocks and other hydrodynamic features are clearly present). Instead, the (thermal) pressure variation is dominated by the mean profile as a function of height above the disk. We briefly examined how well the conditions of hydrostatic pressure equilibrium are fulfilled in these simulations, finding that the thermal pressure is, on average, sufficiently large to support the hot gas (see also Figure 1); however, significant local pressure variations are common, in agreement with previous work (e.g., Joung & Mac Low, 2006).

By the third panel (95 Myr) in Figure 2, the velocity gradient across the kpc-long cloud has caused it to collapse into a relatively narrow high-density region, although with long cloud chains extending down toward the disk. Outflow begins to re-establish by the last time step (105 Myr), and will blow back a small portion of injected gas that did not manage to reach the disk (see also Section 3.4). The pressure, temperature and velocity slices at this point show the first SN shocks behind the inflowing cloud. Movies showing the full frame-by-frame time evolution for all three simulations are provided in the supplementary materials.

Figure 3 shows one-dimensional profiles of the average density, temperature, pressure, and z-velocity as a function of vertical height in each simulation (including a run without inflow). The lighter transparent lines separate the properties of the injected gas (only) from those of the total gas using a normalized colour density cut of $>50\%$. For consistency, these profiles are made when the gas injection is complete at the last time step before the top boundary is switched back to outflow (i.e. the one-kpc cloud has just fully entered the grid). For the BURST50 simulation, this time step corresponds to the first panel (75 Myr) in Figure 2.

In all three simulations, as can be seen in the velocity profiles, the injected gas enters the box and accelerates (gravity only) towards the disk. It quickly decelerates when it reaches $z\sim 1.5-2.0$ kpc, forming a high-density reverse shock. This high-density region is spread between two drops in the velocity profile, beginning with a drop in the injected gas velocity and ending with a drop in the total gas velocity at the boundary where the injected gas and the outflow collide. The vertical extent and amplitude of the density enhancement in the region between the reverse shock and contact discontinuity are inversely related and set by the velocity gradient across the contact discontinuity, which varies between simulations. The 100 km s*-1* (BURST100) run experiences the sharpest deceleration upon impacting the outflowing gas, resulting in the highest density enhancement over a narrow ($\sim 0.25$ kpc) extent. The velocity gradient in the 10 km s*-1* (IMD10) run is much smaller and its broad, low-density reverse shock region spans $\sim 1$ kpc. In all simulations, injected gas that makes its way past the contact discontinuity through low-density gaps in the outflow resumes accelerating towards the disk and fragments into smaller clouds as seen in Figure 2. These clouds are denser and cooler than the average surrounding medium. We discuss these cloud structures and the role of mixing in Section 3.3.

In Figure 4 we show phase diagrams for each simulation at three different time steps. Again, for consistency, we compare across simulations at time steps with similar injected cloud heights. We have superimposed the injected gas properties (colour scale) on top of the total gas properties (gray scale) using the normalized colour density cut. Gas is initially injected in the warm phase with a temperature of $10^{4}$ K and a density of $\rm 8.37\times 10^{-26}\ g\ cm^{-3}$. Over time the inflow distribution broadens to encompass a wide range of densities and temperatures; however, the majority of the injected gas mass cools to lower temperatures (below $10^{4}$ K) and higher densities.

3.3 Cloud Properties

In studying the gas evolution, we find that the injected gas fragments into strands of dense cold cloud structures (Figure 2). In order to study these structures in detail, we create zoomed-in slices of the density, temperature, pressure, z-velocity, colour density, and normalized colour density surrounding the disk. These slices are shown in Figure 5 for the BURST50 simulation. Each slice is centered on the galactic disk with a z-axis extent of $2$ kpc, encompassing the first resolution refinement at $\pm 1$ kpc and the second refinement at $\pm 0.5$ kpc. We have shifted the time step in each panel forward by 10 Myr relative to Figure 2 in order to allow the injected gas time to extend into this region. The cleanest example of these structures is the second time step (95 Myr) in Figure 5, right before the injected gas hits the disk. The density, temperature and velocity slices show that these cloud structures are dense, cold and inflowing towards the disk. We see from the colour density structure that these clouds are composed of injected gas; however, in the normalized colour density structure, it is clear they are not pristine. Recall that the normalized colour density field tracks the fraction of inflowing gas in each cell (so that entirely pristine injected gas has a normalized colour density of 1). The injected cloud structure seen at 95 Myr has a normalized colour density of about [0.6 - 0.8], meaning that the injected gas in these cloud structures has mixed with the outflow/halo gas. The gas surrounding these dense cloud structures has a slightly lower normalized colour density but has also mixed with some of the injected gas. This mixing picture is consistent with the late time phase diagrams in Figure 4, where we see a growing component with higher temperatures and lower densities in the injected gas. With a normalized colour density cut of $50\%$ we are certainly selecting some of this injected gas mixed with the hot, low density outflow.

In Figure 6 we make a quantitative measurement of the mixing between the injected and outflowing gas across simulations. We have applied a z-position cut from z = -0.2 to z = -1 to select gas in the first refinement region above the disk (recall that negative z-position values are above the disk). An additional density cut ($\rm\rho>5\times 10^{-25}\ g\ cm^{-3}$) is applied to each simulation in order to select for gas in the observed cloud structures. Using this cut data, we plot the total mass distribution binned by purity fraction (normalized colour density) and z-velocity. Note that we have not made a purity cut to select only the injected gas; therefore, this plot includes the mass in high density outflow structures as well. In the BURST50 and BURST100 simulations, we find that the majority of the sampled gas lies in two components: a very low-purity outflowing component (with z-velocity $<$ 0) and a higher-purity inflowing, mixed component. Most of this higher-purity component is not pristine, and extends down to a purity fraction of $50\%$ in the BURST50 simulation. We find slightly less mixing in the BURST100 simulation, where the majority of the gas mass does not mix below $\sim 60\%$. In the IMD10 run we only see one component. Both the low-purity outflowing gas and the pristine gas are completely absent and all of the selected gas is mixed with a purity fraction of $\sim 40-90\%$. Note the strong correlation between inflow velocity and purity fraction, an indication that deceleration of the clumps is driven primarily by mixing with the outflowing gas. See also the lower panel in Figure 1; the pressure of the outflowing gas cannot contribute much to the deceleration of the inflow, suggesting an alternate source such as mixing is responsible.

3.4 Inflow Fate

Injected gas that makes it to the galactic disk will either successfully accrete onto the disk or pass through it entirely. We study the fate of the injected gas by splitting the simulation box into three segments. Each segment is a smaller box with the same x-y dimensions as the simulation box and a non-overlapping vertical extent. The middle segment is centered on the galactic disk, extending $\pm 0.2$ kpc in the z-direction. The other two segments are above and below the disk starting at z = $\pm 0.2$ kpc and extending to the edge of the box at z = $\pm 2.5$ kpc. These three segments combined encompass all of the gas in the simulation. We calculate the mass of injected gas in each segment by integrating the product of the colour density and the cell volume throughout each segment. We do not use the normalized colour density cut, as discussed in previous plots, because this quantity includes some of the outflow mixed with the injected gas and would violate our expectation of injected mass conservation as material is exchanged between segments.

Figure 7 shows the integrated mass of the injected gas (left) and the mass injection rate (right) in each segment as a function of time. Note that the total injected gas mass (dashed line) is a roughly constant value equal to the total mass in Table 1. The only time the injected mass is lost from the system is towards the end of the simulation when it either passes through the disk and exits the bottom z-boundary or gets swept up in the re-established outflow and exits the top z-boundary. At the start of each simulation, a finite gas mass has been injected into the top segment above the disk. As the injected gas moves towards the disk it begins spilling over into the disk segment. Any mass that does not stick in the disk then passes through to the bottom segment below the disk. At the finial time step in Figure 7, the injected gas mass in the disk segment has been accreted. Mass remaining above or below the disk has been pushed back towards the top boundary or passed through the disk respectively. We find in Figure 7 that the fate of the injected gas depends strongly on the injection velocity. In the faster 50 km s*-1* and 100 km s*-1* runs, nearly all of the gas makes it to the disk, but only a small fraction (20$\%$ of the total injected gas mass in the BURST50 simulation and 10$\%$ in the BURST100 simulation) sticks in the disk and actually “accretes”. The majority of the injected gas in these high-velocity runs plunges straight through the disk into the bottom segment. This picture is confirmed in the late-stage time evolution panels in Figure 5, where large clouds of injected gas are seen exiting the disk and outflowing towards the bottom z-boundary. In the 10 km s*-1* run, only about half of the injected gas makes it to the disk, but nearly all of that mass (50$\%$ of the injected gas mass) accretes. The mass that does not make it into the disk segment in the IMD10 run is swept up in the re-established outflow and exits the top z-boundary of the box.

In Figure 8 we quantify the temperature of the injected gas accreting onto the galactic disk. The top panel shows the initial temperature distribution of the injected gas in each simulation at a time step directly before the gas hits the disk (panel 3 in Figure 2). We see a spike in all simulations at the injection temperature of $10^{4}$ K. However, in all three simulations the gas quickly redistributes towards lower temperatures. The bottom three segments show the injected gas temperature evolution as a function of time in the disk segment of each simulation. Here we have split the gas into four temperature regimes: cold gas (T $<10^{3}$ K), warm gas ($\rm 10^{3}\leq T<10^{4}$ K), warm-hot gas ($\rm 10^{4}\leq T<3\times 10^{5}$ K), and hot gas (T $\geq 3\times 10^{5}$ K). While a large fraction of warm gas does pass through the disk in the BURST50 and BURST100 simulations, the large majority of the accreted injected gas present in the disk at the end of all three simulation has cooled below $10^{3}$ K.

4 Observational Predictions

We use this suite of simulations to predict observational signatures of gas accretion. These predictions were made using reference tables generated by Corlies & Schiminovich (2016) with the photoionization code CLOUDY (Ferland et al., 1998). Specifically, we are using the g1q1 run at redshift z = 0 detailed in Corlies & Schiminovich (2016). This model was interpolated over a temperature range spanning $[10^{3}-10^{8}]$ K and a number density range spanning $[10^{-6}-10^{2}]\rm\ cm^{-3}$. Density and temperature values from our simulation suite were referenced with these interpolated CLOUDY tables to model the ion fraction and emissivity used in the following analysis. In Section 4.1 we generate mock column density observations in HI and OVI. We chose these two ions because they have observable counterparts and are representative of low and high-energy ionization species. Simiarly, we create surface brightness maps of H$\alpha$ and OVI emission in Section 4.2. The photoionzation model applied assumes an extragalactic UV background (EUVB) as well as solar metallicity and abundances. This model does not include a source of ionizing photons from the disk. We do not expect ionizing radiation from the disk to significantly affect our OVI predictions, as we find collisional ionization is the dominant source of OVI ionization in our simulations. Ionizing photons from the disk may significatly affect the the HI and H$\alpha$ results; however, accurately accounting for radiative transfer of hydrogen ionizing radiation is a beyond the scope of this paper (e.g., Wood et al., 2010).

4.1 Absorption: Column Densities of HI and OVI

In order to calculate the number density $n_{Xi}$ for an element $X$ in ionization state i (e.g. HII: $X=$ H and $i=2$), we substitute the simulation density ($n_{H}$) and the ion fraction computed by CLOUDY $(n_{Xi}/n_{X})$ into Equation 1. For metal ions, we add a factor for the elemental abundance ($n_{X}/n_{H}$) relative to hydrogen, which is the simulated metallicity scaled by the solar abundance (Corlies & Schiminovich, 2016). We create an approximate metallicity field, assuming that the injected gas has a metallicity of $0.2Z_{\odot}$ (Prochaska et al., 2017). The rest of the gas in the simulation volume is given solar metallicity, including the disk, outflow, and disk-halo gas.

[TABLE]

Once we have the ion number density $n_{Xi}$, the projection of $n_{Xi}$ along the line of sight (the x-axis of the simulation box) is the column density shown in Figure 9. Here we compare the BURST50, BURST100, and IMD10 simulations. In addition, we have included an outflow panel (without any inflowing gas) at the time step immediately preceding the gas injection. The right panels in Figure 9 are 2D edge-on column density maps of HI (top) and OVI (bottom). In the left panels, we have averaged across the width of each projection, creating 1D profiles of column density with respect to z-position. These maps reveal a number of interesting structures in the injected gas, outflowing gas, and galactic disk. We compare these predicted column density features with observations. We restrict our comparison to observations of the Milky Way because the initial conditions in this simulation were modeled after the solar neighborhood. In addition, the proximity of studying our own galaxy affords us a wealth of data.

Regions of high HI column density seen in the top panels of Figure 9 probe cold neutral gas. This gas is especially prominent at three locations: (i) the galactic disk, (ii) behind the shock front, and (iii) in the cold high-density cloud structures above the disk. We discuss each of these features in comparison with observations of atomic hydrogen 21-cm emission in the Milky Way.

In the galactic plane, measurements of the HI distribution probe the structure of the galactic disk and the properties of the interstellar medium. From Figure 9, we find a mean log HI column density of $[21.1-21.5]$ at the disk mid-plane across the simulations. We compare this prediction with maps of the integrated HI 21-cm emission intensity from the Leiden-Argentine-Bonn (LAB) all-sky survey (Kalberla et al., 2005). Using the LAB data cube, the mean log HI column density integrated over the velocity range $\pm 400$ km s*-1* within $\pm 10^{\circ}$ of the Galactic plane is 21.5. This number is consistent with our predictions; however it may be a lower limit on the true Galactic column density because it was derived under the assumption that HI 21 cm emission in the Galactic plane is still optically thin.

Beyond the Galactic plane, a number of HI features are observed at the disk-halo interface. In our simulations we find dense cloud structures within about 1 kpc of the disk. As seen in the outflow panel in Figure 9, some of these cloud structures are associated with the SN feedback. The rest are formed from the fragmented injected gas inflowing towards the disk. In Figure 9, both sets of clouds have HI column densities of a few times $\rm 10^{19}\ cm^{-2}$. Similar neutral hydrogen clouds have been observed at the disk-halo interface of the Milky Way in a number of HI surveys. Tangent point observations of 21 cm emission with the Green Bank Telescope (GBT) discovered $\sim 30$ such cloud structures within 0.5 to 1.5 kpc of the disk with HI column densities in the range [0.4 - 2.7] $\rm\times 10^{19}\ cm^{-2}$ (Lockman, 2002). These clouds are co-rotating with the galactic disk, exhibiting only low-velocity deviations on the order of 15 km s*-1*. Recent advancements in observational survey sensitivity, speed, and resolution have led to the cataloging of thousands of neutral hydrogen clouds surrounding the Milky Way. The GALFA-HI survey conducted at the Arecibo Observatory identified 1245 low to intermediate velocity clouds with HI column densities ranging from $[0.1-2.8]\rm\times 10^{19}\ cm^{-2}$ (Peek et al., 2011; Saul et al., 2012). This sample is divided into cold and warm clouds, with median column densities of $\rm 0.6\times 10^{19}\ cm^{-2}$ and $\rm 0.9\times 10^{19}\ cm^{-2}$ respectively. The cloud column densities in our simulations are certainly consistent with the broad range of low-intermediate velocity cloud properties, but they exceed the mean observational column densities. Its possible that this column density excess is caused by the 2 pc minimum grid resolution in our simulations, which constrains clouds from fragmenting below the resolution limit.

At slightly higher latitudes, our simulations predict significant HI column densities behind the shock front in the pile-up of injected gas. This wall feature is produced by the sudden deceleration of the injected gas in contact with the SN-driven outflow. Gas in this density enhancement is able to cool down to (and below) $10^{4}$ K, increasing the neutral fraction and raising the HI column density to $\sim 10^{19}-10^{20}\rm\ cm^{-2}$, as seen in Figure 9. This HI feature is much less pronounced in the IMD10 simulation due to the smaller velocity gradient across the injected gas cloud (and resulting lower densities). Similar high-density wall/ridge structures have been observed alongside high-velocity clouds interacting with the Galactic halo and ISM. We compare with the Smith Cloud, a high-velocity cloud interacting with the Galactic halo. The Smith Cloud is rather large, with an area of $\sim 3\times 1$ kpc and an HI mass of $\geq 10^{6}\rm\ M_{\odot}$; however, the tip of the Smith Cloud is only 2.9 kpc below the galactic plane, making it a rough analogue to the injected gas in our simulations. GBT 21-cm observations of the Smith Cloud show a decelerated ridge of high-density gas $\rm N_{HI}=2\times 10^{20}\ cm^{2}$ along the edge of the cloud (Lockman, 2016; Lockman et al., 2008). The large HI mass of this ridge associates it with decelerated gas from the Smith Cloud and not swept-up material from the surrounding halo gas (Hill et al., 2013). The Smith Cloud also appears to be fragmenting at this boundary, exhibiting indications of Rayleigh Taylor and Kelvin Helmholtz instabilities as it interacts with gas in the Galactic halo (Lockman, 2015)

In contrast to HI, regions of high OVI column density probe highly ionized, usually hot ($10^{5.5}$ K) gas. We find that collisional ionization is the primary ionization source for the OVI gas in our models (see also Smith et al., 2011). Accordingly, the distribution of OVI in our simulations primarily traces the supernovae driven outflow. We compare our OVI column density predictions with Li et al. (2017b), who performed a detailed study of the OVI density distribution in simulated outflow generated using an identical SN prescription. In Li et al. (2017b), the OVI density was calculated using the product of the O/H solar abundance ($5\times 10^{-4}$) and the temperature-dependent OVI ionization fraction from Sutherland & Dopita (1993). The OVI column densities stated in Li et al. (2017b) were calculated for a time-step of 140 Myr with error bars noting the time variation in the latter half of the simulation. Our outflow time step of 50 Myr is therefore outside the range of these error bars. Although we cannot make a direct comparison between the two sets of the predictions, we note similarities and differences below. Using CLOUDY, we predict an average OVI density of $n_{OVI}=1.8\times 10^{-8}\rm\ cm^{-3}$ within $\pm 100$ pc of the disk. This is roughly consistent (to within $2\sigma$) with the average OVI density of $n_{OVI}=(2.4\pm 0.3)\times 10^{-8}\rm\ cm^{-3}$ predicted by Li et al. (2017b) within $\pm 100$ pc of the disk. When we calculate the average OVI density within $\pm 200$ pc from the disk at 50 Myr we find that the value remains approximately constant at $1.8\times 10^{-8}\rm\ cm^{-3}$. Li et al. (2017b) reports a very similar average OVI density of $n_{OVI}=(1.8\pm 0.3)\times 10^{-8}\rm\ cm^{-3}$ within $\pm 200$ pc at 140 Myr, which has dropped relative to their prediction close to the disk. All of these values are consistent with absorption line studies of nearby stars, which measure the OVI density of the local ISM to be $n_{OVI}=(2.2\pm 1.0)\times 10^{-8}\rm cm^{-3}$ (Jenkins, 1978; Oegerle et al., 2005; Savage & Lehner, 2006; Bowen et al., 2008; Barstow et al., 2010).

4.2 Emission

We next turn to line emission as a complementary probe of the disk-halo conditions, again for the same two ions. The emission intensity ($\Gamma$) is calculated from the CLOUDY emissivity using Equation 2. Here $\epsilon(\rho,T)$ is the temperature and density dependent emissivity, $h$ is Planck’s constant, $c$ is the speed of light, and $\lambda$ is the wavelength of the recombination transition ($\lambda 6563\AA$ for $H\alpha$ and $\lambda 1032\AA$ for OVI). The projection of the emission intensity integrated along the line of sight is the surface brightness shown in Figure 10. On the right we have plotted edge-on projections of the $\rm H\alpha$ and OVI emission surface brightness for each of the three simulations and compared them with surface brightness predictions for an outflow time step preceding the gas injection. The corresponding average emission profiles as a function of z-position are plotted on the left panels.

[TABLE]

The H$\alpha$ emission in Figure 10 traces high-density structures of warm, partially ionized gas and highlights regions of strong gas cooling. Most notably, regions of high H$\alpha$ emission intensity correspond to (i) the high-density ISM ejected from the disk via supernovae feedback, and (ii) the dense wall of injected gas piling up behind the shock front.

We compare these predictions with observations from the Wisconsin H$\alpha$ Mapper (WHAM) (Haffner et al., 2003). WHAM surveyed the distribution and kinematics of H$\alpha$ emission in and around the Milky Way. These observations probe the warm/warm-hot gas in the ISM as well as revealing a number of H$\alpha$ structures at the disk-halo interface. In Figure 10 the galactic disk and supernova-driven outflows produce H$\alpha$ emission intensities ranging from $\sim 10^{3}$ to $10^{5}\rm\ photons\ s^{-1}\ cm^{-2}\ sr^{-1}$. Gaps in the inhomogeneous outflow reveal the hot halo gas with essentially no H$\alpha$ emission intensities. As with the column density, we expect to under-predict the emission intensity due to the shortened projection path length (0.35 kpc) through our simulation box. WHAM observations along the Sagitarius-Carina Arm reveal similar H$\alpha$ emission values ranging from $>2\times 10^{6}\rm\ photons\ s^{-1}\ cm^{-2}\ sr^{-1}$ in the galactic plane to $\sim 1\times 10^{4}\rm\ photons\ s^{-1}\ cm^{-2}\ sr^{-1}$ beyond the disk (Krishnarao et al., 2017). Models of the disk scale-height place these emission measurements within a vertical height of 3 kpc of the Galactic disk, making them an excellent analog for our simulation box height. The WHAM emission maps also reveal patches below the detection limits similar to the predicted gaps in our H$\alpha$ emission maps.

Another prominent H$\alpha$ feature in our simulation suite is the dense wall of decelerated injected gas. In Figure 10 we see a thick wall of elevated H$\alpha$ emission trailing the leading edge of the injected gas. Looking closely, this wall is actually composed of many thin density enhancements projected on top of one another in the edge-on inclination. This effect is due to the inhomogeneous outflow halting the injected gas at different heights above the galactic disk. As mentioned in Section 4.1, we compare this feature with the Smith Cloud, an HVC at the Milky Way disk-halo interface interacting with the ISM. More than half the mass of the Smith Cloud is ionized, producing a smooth distribution of H$\alpha$ emission (Wakker et al., 2008; Hill et al., 2009). In addition, H$\alpha$ observations with WHAM reveal a ridge of elevated H$\alpha$ emission trailing the HI ridge by a few hundred pc. Although the HI and H$\alpha$ ridges are not coincident, their similar morphology suggests a physical connection. Gas in this ionized feature is likely shock heated, with an H$\alpha$ emission intensity ranging from $[3.17-5.56]\times 10^{4}\rm\ ph\ s^{-1}\ cm^{-2}\ sr^{-1}$ (Hill et al., 2013). In our simulations, H$\alpha$ emission associated with the deceleration of injected gas ranges in mean intensity from about $[1-15]\times 10^{4}\rm\ ph\ s^{-1}\ cm^{-2}\ sr^{-1}$ (Figure 10). Overlapping H$\alpha$ ridges are seen trailing the cloud front by as much as half a kpc. While we cannot make a direct comparison with the neutral and ionized ridge morphology observed in the Smith Cloud, we note that we are able to reproduce a reasonable range of elevated H$\alpha$ emission intensities where the shocked ionized component is separated from the decelerated neutral gas.

Similar to H$\alpha$, OVI emission also traces gas cooling, but in the highly ionized hot gas. In agreement with the low OVI column densities predicted in the injected gas in Figure 9, we do not predict significant OVI emission from the injected gas in Figure 10. We conclude that the injected gas is not heated (above $10^{5.5}$ K) and re-cooled on its path the the disk. However, we do predict significant OVI emission from the SN-driven outflow with an average intensity of $1.84\times 10^{4}\rm\ ph\ s^{-1}\ cm^{-2}\ sr^{-1}$ projected along the z-axis. We compare our Cloudy prediction with OVI emission predictions from Li et al. (2017b) which were calculated using the OVI ion fraction described in Section 4.1 and an oscillator strength of 5 for the $\lambda 1038\AA$ emission line (Osterbrock, 1989). As noted in Section 4.1, Li et al. (2017a) makes emission predictions at a time step of 140 Myr with error bars time averaged over the later half of the simulation. Our final outflow time step at 50 Myr falls outside the range of these error bars; however, we still make an approximate comparison. The z-axis projection in Li et al. (2017b) predicts an OVI $\lambda 1032\AA$ emission intensity of $[0.5-6]\times 10^{3}\rm\ ph\ s^{-1}\ cm^{-2}\ sr^{-1}$ at a time step of 140 Myr, which is about a factor of three lower than our emission predictions. Note that both our predictions and the Li et al. (2017b) predictions are the result of projecting through the entire simulation box (along the z-axis), and need to be divided by two in order to compare with the following observations taken within the Milky Way. Observations with FUSE detected OVI $\lambda 1032\AA$ emission along 23 sight lines with intensities ranging from $(1.9-8.6)\times 10^{3}\rm\ ph\ s^{-1}\ cm^{-2}\ sr^{-1}$ (Otte & Dixon, 2006). Despite extinction, higher emission intensities in this survey were actually found at low-latitude sight lines in the nearby ISM. In addition, high-latitude observations of OVI emission have been obtained with SPEAR close to the North Galactic Pole, measuring an emission intensity of $(0.5-2.0)\times 10^{4}\rm\ ph\ s^{-1}\ cm^{-2}\ sr^{-1}$ (Welsh et al., 2007). Our OVI emission predictions (as well as the prediction from Li et al. (2017b)) are roughly consistent with this range of observational values.

5 Discussion

In this work, our inflowing gas cloud was deliberately idealized: we assumed a constant density, temperature and inflow velocity. This was done in part because our goal is to understand the physical processes operating in the inflow/outflow interaction, and our outflow properties are already quite complicated, with multiphase gas ejected from realistically generated supernovae explosions. We briefly discuss physical origins for this inflowing gas cloud in Section 5.1 and elaborate on select results from Section 3 including the evolution of the injected cloud structures (Section 5.2) and the flow of injected gas plunging through the galactic disk (Section 5.3). We compare this work with the existing literature (Section 5.5) and make note of some of the simulation limitations and assumptions (Section 5.6)

5.1 Inflow Origin

The original source of the injected gas cloud in this simulation could be either cosmological inflow, gas stripping from satellites, or gas recycling, all of which are expected to fuel the galactic disk. Cosmological inflow in the form of substantial cold filamentary flows has been seen in a wide variety of simulations (e.g., Kereš et al., 2009; Nelson et al., 2013; Hobbs et al., 2013; Joung et al., 2012b). The properties of these flows close to the disk-halo interface are not entirely clear. As demonstrated by the breadth of work discussed in Section 5.5, understanding the flow of cold clouds through a hot CGM is challenging and requires additional numerical, observational and analytic studies. However, it seems possible and maybe even likely that at least some gas clouds with a cosmological origin can get close to the disk without being completely mixed into the CGM (Stewart et al., 2017).

Additionally, the inflow in these simulations could originate from recycled gas flow. Assuming that outflows are highly mass loaded (as seem to be required to reproduce galaxy properties), much of the ejected gas will fall back into the CGM, and eventually back onto the galactic disk. This was found in early (relatively low-resolution) simulations of momentum-driven winds (Oppenheimer et al., 2010) but has since been confirmed and quantified in higher resolution simulations with more physically based feedback (Muratov et al., 2015; Anglés-Alcázar et al., 2017).

5.2 Cloud Evolution

Despite our simplified initial conditions, we found that the injected gas rapidly fragmented into small clumps through a series of instabilities, both hydrodynamic (Rayleigh-Taylor, Kelvin-Helmholtz) and thermal (via radiative cooling). These clumps initially accelerated ballistically from their starting inflow velocity (see Figure 3). The amount of acceleration immediately after injection (before slowing due to interactions with the outflowing material) was larger for simulation runs with a lower injection velocity (e.g. the IMD10 gas gained 40 km s*-1* while BURST100 only added 20 km s*-1* to its initial velocity).

Over time, the inflowing clouds in all three simulation runs were decelerated as they interacted with the outflowing gas. As demonstrated in Figure 6, clouds within 1 kpc of the disk have a narrow velocity distribution that appears to peak around a typical cloud velocity in each simulation ($\sim$ 40 $\rm km\ s^{-1}$ in the IMD10 simulation to $\sim$ 100 km s*-1* in the BURST100 simulation). It is tempting to associate this with a “terminal” velocity; however, the analogy with the classic computation should not be taken too far, as the clouds are not coherent entities. Unlike a penny dropped from a tall building, the clouds have no interstitial forces to hold them together. Instead, as evident from the tight correlation in Figure 6, deceleration is directly associated with gas mixing. Therefore, a simple momentum balance argument is probably a better explanation of the resulting cloud velocity distribution.

Irrespective of the physics causing the cloud deceleration, it appears that the fragmentation of our injected gas into dense ($n\sim 1$ cm*-3*), cold ($T<10^{4}$ K) clumps with velocities of 50-100 km s*-1* is difficult to avoid. Even in our IMD10 case, where the inflow ram pressure is below the outflow ram pressure, we find inflowing clouds penetrating down to the disk with velocities of  40 km s*-1*. In this low-velocity case, fragmentation of the injected cloud results in dense cloud clumps that are able to overcome the outflowing ram pressure. Radiative cooling is key to the creation of these cold dense clouds and should be traced by the associated Ly$\alpha$ and H$\alpha$ emission (as described in Section 4).

Another interesting aspect of the cloud evolution is the spatial correlations in the clouds. This can be seen in Figure 5 but is much clearer in the supplied movies. The clumps are not randomly distributed throughout the hotter background material filling the disk-halo interface, but tend to occur in chains or filaments (which themselves are clustered). This clustering appears to be due to a combination of effects. The physical process of cloud fragmentation causes larger clouds to break up into smaller localized groups, and this group morphology facilitates a hydrodynamic shielding effect (like Canadian geese flying in formation) that prevents stripping and disruption.

5.3 Outflow Past the Disk

As mentioned in Section 3.4, some of the injected gas passes through the galactic disk and continues towards the bottom boundary of the simulation box. This is especially prevalent in the BURST50 and BURST100 simulations where the inflowing gas has higher ram pressure and is able to penetrate through low surface density regions in the disk. The bottom panels in Figure 5 show the inflow colour density and normalized colour density for gas within $\pm 1$ kpc from the galactic disk in the BURST50 simulation. Recall that the inflow colour density is the density of the injected gas and the normalized colour density is the ratio of injected gas density to total gas density. The third (105 Myr) and fourth (110 Myr) inflow colour density time steps in Figure 5 show the 50 km s*-1* injected gas plunging through the galactic disk. This gas is then swept away from the galactic plane with the SN outflow. Before hitting the disk, the injected gas has a high normalized colour density ($>0.6$). Some mixing with the outflowing gas has occurred, but the injected gas is still close to pristine. After exiting the disk, the injected gas exhibits a broader range of normalized colour density values. We find clumps of injected gas with high normalized colour density values which have not mixed significantly during their passage through the disk. However, we also find injected gas with low ($<0.6$) normalized inflow colour density values, suggesting that some of the injected gas passing through the galactic disk has mixed with the ISM.

5.4 Inflow Extent

In the simulations above, inflowing gas was injected across the entire cross-sectional area of the top z-boundary. This input geometry is somewhat arbitrary, and we have completed an additional simulation with a substantially smaller cloud to ensure that the extent of the cloud does not significantly influence the physical results. This small cloud simulation (SMALL50) uses the same density, temperature, and injection velocity as the BURST50 run. We inject gas from a smaller square region in the center of the top z-boundary amounting to a quarter of the total box cross-sectional area. The rest of the area on the outer edges of the top z-boundary remains set to outflow. For consistency with the other runs, we inject this cloud until it reaches a kpc in length at which time we revert the full cross-sectional area of the top z-boundary to outflow. Due to the reduced cloud extent, the resulting injected cloud mass is only about a quarter of the BURST50 injected cloud mass. The time evolution of the SMALL50 simulation in density, z-velocity, and normalized colour density is shown in Figure 11. We run this simulation to 90 Myrs, by which time the cloud is clearly on its way to interact with the disk. It has fragmented into smaller clumps and descended to within a kpc of the disk midplane.

The structure and evolution of this smaller injected cloud is similar to that of the injected gas in the BURST50 simulation. The injected gas forms a high-density front in contact with the outflowing gas and fragments into dense cold clumps. The maximum z-position of the injected gas in the SMALL50 simulation slightly lags the BURST50 injected gas at each time step. Injected gas in the BURST50 simulation was able to penetrate down towards the disk through gaps in the outflow along the edge of the simulation box where gas has not been injected in the SMALL50 run. The phase diagrams for the SMALL50 simulation have an analogous distribution to the phase diagram of the BURST50 simulation. As in the BURST50 simulation, the evolution of the SMALL50 phase diagrams show the majority of the injected gas redistributing to lower temperatures and higher densities over time. At equivalent time steps, the SMALL50 injected gas is not as close to the disk as the BURST50 simulation. Accordingly, the phase diagram of the SMALL50 injected gas has not fully populated the higher pressure isobars that BURST50 injected gas achieves at lower latitudes. However, this is largely a timing issue (the SMALL50 gas infall occurs more slowly because it was injected directly into a stream of outflowing gas at the box center) and as time goes on, it evolves in the same manner as the BURST50 run.

One notable difference in this geometry is that the outflowing gas is able to pass around the injected cloud creating shearing between the edge of the injected cloud and the adjacent outflowing gas. This shearing strips away a small fraction of the injected cloud mass. Most of this stripped material sits at high-latitudes stalled behind the injected cloud, but some ($<4\%$ of the total cloud mass) is pushed out of the box with the outflowing gas. The normalized colour density field (Figure 11) shows mixing along the shear layer, with the high-latitude stripped gas exhibiting a purity of less than $60\%$.

5.5 Comparison to previous work

We are not the first to explore the evolution and fate of gas flowing into a galaxy. Heitsch & Putman (2009) examined the evolution of similar gas inflow both in simplified wind-tunnel experiments and in hydrostatic stratified atmospheres without outflows. They found that injected gas clouds rapidly generated head-tail morphologies and ablated, generally surviving only a few kpc. This result was dependent on the size of the cloud, with larger clouds exhibiting longer survival times over greater distances. They note the importance of cooling, as some of the warm ionized material stripped from the infalling clouds was later able to condense and re-cool close to the disk. Our results are not inconsistent with this work. The closest comparison is with their larger mass, stratified cases where they find similar cloud fragmentation and cooling due the compression of these cloud fragments on their way through the halo. This cooling increases the HI mass and helps stabilize the clouds against disruption, allowing them to travel greater distances. In agreement with our cloud velocity values, they also showed an early phase of near-ballistic acceleration followed by deceleration when they reached velocities of $\sim 100\rm\ km\ s^{-1}$.

Similar results were also found in a set of simulations carried out using the FLASH code (Kwak et al., 2011; Shelton et al., 2012; Gritton et al., 2014, 2017), where they studied individual cold clouds interacting with a hot atmosphere. Outflows in these simulations were generally implemented using the constant density and velocity wind-tunnel model. These studies also concluded that low mass clouds were quickly destroyed, with larger clouds exhibiting longer survival times. Additionally, they note the importance of cooling and mixing in models of massive clouds ($10^{5}-10^{8}\rm M_{\odot}$), which increase their HI mass on their way to the disk by accreting a substantial amount of ambient halo gas through condensation.

In addition, Marinacci et al. (2010) and Marinacci et al. (2011) examined the evolution of clouds moving through the hot halo near the disk under similar conditions to those examined here. However, in these simulations the clouds were presumed to be fountain flows, the motion was not primarily vertical, and the hot halo was uniform and static. They carried out two-dimensional simulations, likewise reporting that the clouds compress, fragment, and cool efficiency thorough radiative cooling. They also concluded that cold clouds interacting with a hot halo medium will gain mass from their surroundings and decelerate due to the transfer of momentum that accompanies this exchange of material. Thermal conduction was included in the Marinacci et al. (2011) simulations byArmillotta et al. (2016), who found that the conduction in higher temperature halos (above $4\times 10^{6}$ K) suppressed hydrodynamic and thermal instabilities and reduced the importance of re-cooling.

5.6 Simulation Caveats

Inflow conditions: As noted earlier, our inflow conditions are quite simplified. The constant density and velocity of material injected across the upper boundary condition of the simulation box corresponds to large, diffuse clouds. In addition, we have assumed they are relatively cold. Although these bulk conditions are motivated by cosmological simulations, their smooth injection is probably unrealistic. More likely, the clouds would have internal structure related to velocity, density and temperature perturbations (including turbulence). We argue that this would probably only increase the likelihood for the clouds to fragment into smaller clumps, and so strengthen the results that we find here. However, we recognize that very strong turbulent velocities or greatly enhanced temperatures may provide a source of non-thermal pressure and change the cloud evolution entirely.

Radiative transfer: We make quite simplified assumptions about the radiative conditions in the cloud. In particular, we assume optically thin radiative cooling (likely reasonable given the densities) and a uniform photo-heating rate discussed in more detail in Li et al. (2017a). In reality, some ionizing radiation probably leaks out of the disk and may heat the gas (e.g., Bland-Hawthorn & Maloney, 1999; Dove et al., 2000). Also, radiation pressure from both ionizing and infrared radiation could provide additional lift for the clouds (Murray et al., 2010). We suspect both of these effects will be relatively minor in determining the dynamics of the clouds in this simulation. Although heating of low temperature clouds may be important in some cases, radiation pressure is almost certainly sub-dominant to thermal and ram pressure.

Magnetic fields: There is only limited information on the magnetic field strength in clouds above the disk; however, fields are certainly observed to have strengths similar to the thermal pressure in the interstellar medium (e.g., Heiles & Crutcher, 2005). Numerical models of the disk which do include fields (e.g., de Avillez & Breitschwerdt, 2005; Hill et al., 2012; Walch et al., 2015) generally find their dynamical impact to be relatively mild, although there is some indication they may provide additional support for the vertical scale height. If important in inflowing clouds, they may help to prevent fragmentation by providing dynamical support, as well as preventing or slowing the development of the Kelvin-Helmholtz or Rayleigh-Taylor instabilities.

Numerical Resolution: The numerical resolution in our simulations (2/4/8 pc within 0.5/1.0/2.5 kpc of the mid-plane) was chosen to be sufficiently high to resolve the Sedov phase of our SNe explosions within all but the highest density regions Li et al. (2017a); Joung et al. (2012b). However, the cold clouds in the multiphase medium are not fully resolved. Indeed, without the addition of thermal conduction it is unlikely that we can fully resolve the cloud properties at any resolution. Instead they will continue to break up into smaller clumps down to the grid scale (Choi & Stone, 2012).

6 Summary

We have carried out a set of numerical simulations modeling a cloud of gas inflowing towards the galactic disk from a height of 2.5 kpc. We include realistic heating and outflows in order to better understand gas interactions at the disk-halo interface. Our galactic model is taken from Li et al. (2017a) and is aimed at approximating conditions similar to the local neighborhood of the Milky Way. We include SN explosions resolved with a 2 pc spatial grid, which should be sufficient to model their energetic input. The gas accretion is simulated with a smooth inflow of gas at the upper boundary condition. We explore three models which vary in the velocity of the inflow, ranging from 10 to 100 $\rm km\ s^{-1}$. In each case, a similar amount of total gas is injected before the inflow is curtailed and the resulting flow is evolved for  50-100 Myr. We summarize the following key results:

Cloud Evolution The initial diffuse cloud interacts with the outflowing gas, forming a forward shock, contact discontinuity and reverse shock. Most of the cloud mass ends up behind the discontinuity, cooling and fragmenting via Rayleigh Taylor and Kelvin Helmholtz instabilities. The clumps have typical densities of 1 cm*-3* and temperatures of about 1000 K. The dense clumps can accrete onto the disk even when the ram pressure of the outflow exceeds that of the original diffuse cloud (Figures 1, 2 and 4). 2. 2.

Cloud Velocity and Mixing The clouds initially accelerate ballistically before being decelerated by interactions with the outflowing material. Close to the disk, we find typical cloud velocities ranging from 40-100 $\rm km\ s^{-1}$, depending on their initial injection rate. This approximately corresponds to intermediate velocity clouds. The deceleration appears to be largely controlled by mixing, with velocities of dense gas correlating strongly with the amount of mixing the cloud has experienced. Clouds do mix somewhat with the ambient medium but retain more than 50% of their original gas (Figures 3 and 4). 3. 3.

Cloud Fate and Accretion In these models, infalling gas clouds often penetrate through the disk, either because they hit an area with low disk surface density, or because they have a high ram pressure. Most of the penetrating gas is reaccelerated by the outflow on the other side of the disk and expelled. Only in the case with the lowest initial velocity (10 km/s) does a majority of the cloud mass end up accreting on to the disk (Figure 7). 4. 4.

Observational comparison We generate mock observables for HI (a low ionization ion) and OVI (a high ionization ion). We predict that inflow produces significant HI column densities of $\sim 10^{20}$ cm*-2* in the compressed injected gas, with even higher values for simulations with higher inflow rates. This is consistent with observational column densities of IVC’s as discussed in section 4.1. We find correspondingly high H$\alpha$ emission in the high-density (reverse shock) region with an elevated emission intensity of $[1-15]\times 10^{4}\rm\ ph\ s^{-1}\ cm^{-2}\ sr^{-1}$. Closer to the disk, emission from the fragmented cloud structures is difficult to distinguish from the background outflow emission intensity. Because the clouds stay largely cool, we find that there is relatively little OVI absorption or emission. We make comparisons to observations, finding rough agreement with many existing Milky Way measurements. (Figures 9 and 10).

This work highlights the importance of interactions between inflowing gas clouds and outflowing gas at the disk-halo interface. The compression, fragmentation and cooling of inflowing gas results in the formation of cold, dense clouds which descend to the disk (and beyond) without being destroyed. Future work is required to explore a wider parameter space, both of cloud properties and of disk conditions. In addition, we have neglected a number of physical effects, the importance of which should be examined in later work.

We thank Mary Putman, David Schiminovich, Yong Zheng and our referee for their useful discussions and thoughtful feedback, which have helped us improve various drafts of this manuscript. Many thanks to Lauren Corlies for sharing her Cloudy models used in our mock observations. We acknowledge support from NSF grants AST-1312888, AST-1615955 and DGE-1644869, and NASA grant NNX15AB20G, as well as computational resources from NSF XSEDE, and Columbia University’s compute cluster. Simulations were run on the Texas Advanced Computing Center Stampede supercomputer and were performed using the publicly-available Enzo code (http://enzo-project.org) and analyzed with the yt package (Turk et al., 2011). Enzo is the product of a collaborative effort of many independent scientists from numerous institutions around the world. Their commitment to open science has helped make this work possible. The Flatiron Institute is supported by the generosity of the Simons Foundation.

Appendix

Movies for the BURST50, BURST100, and IMD10 simulations can be found at the following url: https://doi.org/10.7916/d8-vgnn-7n73. They are also embedded in the Appendix of the ApJ publication of this work. We have included the corresponding still frames (the first frame of each movie) below. In previous figures we changed the aspect ratio of the simulation slices (compressing them along the z-direction) to ensure they fit in the figure. In these movies/slices we maintain the correct aspect ratio (10:1).