Introducing CGOLS: The Cholla Galactic Outflow Simulation Suite

Galactic outflows are recognized as an important driver of galaxy evolution. Outflows regulate the mass of galaxies and their baryon to dark matter ratios and enrich the circumgalactic and intergalactic media (CGM and IGM) with metals (Dekel & Silk 1986; Katz et al. 1996; Springel & Hernquist 2003; Mashchenko et al. 2008; Oppenheimer & Davé 2008; Scannapieco et al. 2008; Davé et al. 2011; Hopkins et al. 2014; Vogelsberger et al. 2014; Schaye et al. 2015). The input of mass and energy from massive stars, particularly via supernovae, is thought to be the major driver of galactic outflows in low-mass star-forming galaxies (Heckman et al. 1990; Hopkins et al. 2012). Understanding the process by which stellar feedback drives gas out of galaxies is therefore a fundamental goal of current galaxy evolution theories.

Numerical simulations have become a critical tool in our theoretical understanding of galaxy evolution, as they can model and predict nonlinear processes that cannot be easily investigated analytically. Supernova-driven galactic winds are certainly such a process, given their complex observed structures and the number of relevant physical processes governing their evolution. Hydrodynamics, radiation, magnetic fields, conduction, and gravity may all play important roles, and the effects of many of these processes have been investigated via 2D simulations (e.g., Mac Low et al. 1989; Strickland & Stevens 2000), 3D simulations of the inner region of starburst galaxies (e.g., Cooper et al. 2008; Sarkar et al. 2015), global galaxy simulations with static mesh refinement (e.g., Fielding et al. 2017a), local box simulations of dense clouds within a rarefied hot wind (e.g., Schiano et al. 1995; Marcolini et al. 2005; Melioli et al. 2005; Cooper et al. 2009; McCourt et al. 2015; Scannapieco & Brüggen 2015; Brüggen & Scannapieco 2016; Schneider & Robertson 2017), and simulations of patches of a galaxy's interstellar medium (ISM) and regions above the disk (e.g., Joung & Mac Low 2006; Kim et al. 2013; Walch et al. 2015; Martizzi et al. 2016; Li et al. 2017).

Simulating starburst-driven winds has its own challenges, however. In addition to the computational expense of including the relevant physics, the range of spatial scales involved is daunting. Critical parameters of winds include their mass and energy loading, $\dot{M}={\dot{M}}_{\mathrm{wind}}/{\dot{M}}_{\mathrm{SFR}}$ and $\dot{E}=\dot{E}\mathrm{wind}/{\dot{E}}_{\mathrm{SN}}$, which are set on the scale of individual supernova bubbles on the order of a few parsecs (Kim & Ostriker 2014; Martizzi et al. 2015; Kim et al. 2017). The global distribution of sources throughout the disk may have a large effect on the subsequent evolution of the outflows and requires global galaxy simulations to capture (Martizzi et al. 2016; Vijayan et al. 2018). Following the winds out into the CGM makes the problem even more challenging, as the scales now approach the virial radius of a galaxy (∼100 kpc), while individual cool clouds may still exist with size scales ≪1 kpc (Fielding et al. 2017b; McCourt et al. 2018).

While the challenges in dynamic range are substantial, scientific access to ever-increasing computing power continues to push the boundaries of possibility. In this paper, we describe the Cholla Galactic OutfLow Simulations (CGOLS) project, a suite of simulations that attempts to tackle the problem of understanding galactic winds via a systematic approach. Beginning with a simple, analytically tractable feedback model for a galactic wind, we progressively increase the physical complexity of the feedback implementation and the included physics through each simulation in the suite, allowing us to investigate for each case the way in which the numerical results diverge from simpler theoretical expectations. In addition to a novel setup, a primary advancement the CGOLS suite contributes to the study of these winds is the extreme resolution made possible through the use of our GPU-based simulation code, Cholla (Schneider & Robertson 2015). While we cannot yet achieve 1 pc resolution over a volume of 100 kpc³, the production simulations described in this paper each have a resolution of ≈5 pc over a volume of 10 × 10 × 20 kpc³, which represents over an order-of-magnitude improvement as compared to any other isolated-galaxy simulations of these phenomena (see Figure 1). This vast resolving power allows us to tackle the problem of simulating galactic winds with increased fidelity and enables us to pay careful attention to the role that finite numerical resolution may play in our understanding of galactic winds via simulation.

Zoom In Zoom Out Reset image size

Broadly, CGOLS consists of a set of global galaxy simulations designed to study the evolution of galactic outflows between radii of 1–10 kpc. The GPU-native nature of Cholla allows us to run on Titan, as of 2017 November, the fastest computer available for public science in the United States and the fifth fastest in the world.⁴ The immense power of Titan and the efficiency of Cholla make it straightforward to achieve resolutions for these isolated-galaxy simulations that far outstrip those previously possible. In order to place these simulations in context, we plot in Figure 1 a comparison between the total hydrodynamic volume elements (particles or cells) and resolution of the CGOLS suite and those of other simulations that have been used to study the structure and dynamics of galactic outflows. The relative locations of a representative sample of ISM patch and isolated-galaxy simulations, zoom-in and static mesh refinement simulations, and cosmological simulations are outlined. As this plot demonstrates, the CGOLS project investigates winds in an entirely new region of parameter space, one in which the effects of global wind geometry and propagation into the CGM can be explored, while the resolution remains sufficient to better capture small-scale hydrodynamic and/or thermal instabilities that affect the multiphase structure of the outflow. Essentially, we are applying the numerical resolution of the current premier cosmological simulations to the inner halo of a single galaxy.

We organize this paper as follows. In Section 2, we provide details of the simulation suite, including the setup of the isolated galaxy and halo, the feedback model, and the relevant details of our code. In Section 3, we present results from the first simulation in the CGOLS suite, which investigates an adiabatic, nuclear starburst wind model. In Section 4, we connect this simulation to observations of the nuclear starburst galaxy, M82, with a focus on the observed diffuse soft X-ray emission. In Section 5, we explore the role of resolution and discuss convergence of the results. We conclude in Section 6. A companion paper (Schneider et al. 2018) presents results from simulations in the suite that include radiative cooling, and therefore investigates the potential multiphase nature of these outflows.

Table 1 lists the base set of simulations in the CGOLS suite. The scientific results presented in this paper primarily result from the first three simulations shown in the table, but Table 1 includes additional CGOLS simulations for completeness. Each CGOLS simulation shares the same galaxy initial conditions and halo model, as described below.

2.1. Simulation Setup and Galaxy Initial Conditions

A schematic representation of the simulation volume is shown in Figure 2. Centered in the volume is a rotating gas disk, initially surrounded by a static, adiabatic halo. The first feedback model we study in this project is simple and consists of a constant injection of mass and thermal energy into the central region of the galaxy, which drives a collimated wind. The fiducial simulation volume spans a region 10 × 10 × 20 kpc³, with 2048 × 2048 × 4096 cells, giving a constant resolution of Δx = 4.9 pc across the entire volume. Additional simulations are also performed at lower resolution (see Table 1).

Zoom In Zoom Out Reset image size

2.1.1. Gravitational Potential

We use a static gravitational potential to model the contributions of both the stellar component of the galactic disk and dark matter halo throughout the simulations. Both the disk and halo are modeled after the gas-rich nearby starburst galaxy, M82. Self-gravity of the gas is not included in these initial CGOLS simulations. We use a Miyamoto-Nagai profile for the galaxy's stellar disk component (Miyamoto & Nagai 1975), defined by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{stars}}(r,z)=-\displaystyle \frac{{{GM}}_{\mathrm{stars}}}{\sqrt{{r}^{2}+{({R}_{\mathrm{stars}}+\sqrt{{z}^{2}+{z}_{\mathrm{stars}}^{2}})}^{2}}},\end{eqnarray} \tag{ 1 }$$

where r and z are the radial and vertical cylindrical coordinates, M_stars = 10¹⁰ M_⊙ is the mass of the stellar disk (Greco et al. 2012), R_stars = 0.8 kpc is the stellar scale radius (Mayya & Carrasco 2009), and z_stars = 0. 15 kpc is the stellar scale height (Lim et al. 2013). We use an NFW profile to represent the halo potential (Navarro et al. 1996), defined by

$$\begin{eqnarray}&&{{\rm{\Phi }}}_{\mathrm{halo}}(r)=-\displaystyle \frac{{{GM}}_{\mathrm{halo}}}{r[\mathrm{ln}(1+c)-c/(1+c)]}\mathrm{ln}\left(1+\displaystyle \frac{r}{{R}_{\mathrm{halo}}}\right),\end{eqnarray} \tag{ 2 }$$

where r is the radius in spherical coordinates, M_vir = 5 × 10¹⁰ M_⊙ is the dark matter mass of the halo, c = 10 is the halo concentration, and R_halo is the scale radius of the halo, which is set by R_h = R_vir/c = 5.3 kpc with our estimate of the virial radius, R_vir = 53 kpc. We have estimated the dark matter mass and virial radius of the halo based on scaling relations between the stellar content of galaxies and their dark matter halos (e.g., Kravtsov 2013), but expect values within a factor of a few for the halo parameters to have little effect on the results of our simulations. The total gravitational potential in our simulations is then Φ = Φ_stars + Φ_halo. Figure 3 shows the theoretical stellar rotation curve in the x–y plane of our simulated galaxy potential plotted against stellar data from M82 (Greco et al. 2012). (The initial simulations in the CGOLS suite do not include star particles and only model gas. Correspondingly, Figure 3 merely demonstrates that the adopted potential provides a good match to the M82 observations.)

Zoom In Zoom Out Reset image size

2.1.2. Disk Gas

The galaxy itself in our simulations consists solely of gas. The gas is distributed radially with an exponential surface density profile defined by the scale radius, R_gas = 1.6 kpc = 2 × R_stars, and total gas mass, M_gas = 2.5 × 10⁹ M_⊙ (Greco et al. 2012), according to

$$\begin{eqnarray}&&{\rm{\Sigma }}(r)={{\rm{\Sigma }}}_{0}\exp \left(-\displaystyle \frac{r}{{R}_{{\rm{g}}}}\right),\end{eqnarray} \tag{ 3 }$$

where r is the cylindrical radial coordinate, and Σ₀ is the central surface density, ${{\rm{\Sigma }}}_{0}={M}_{\mathrm{gas}}/(2\pi {R}_{\mathrm{gas}}^{2})$. The disk gas is isothermal with temperature T_disk = 10⁴ K and is set in vertical hydrostatic equilibrium with the gravitational potential. This is achieved by solving for the vertical density distribution at each radial location,

$$\begin{eqnarray}&&\rho (z)={\rho }_{0,{\rm{d}}}\exp \left[-\displaystyle \frac{{\rm{\Phi }}(z)-{{\rm{\Phi }}}_{0,{\rm{d}}}}{{c}_{{\rm{s}},{\rm{d}}}^{2}}\right],\end{eqnarray} \tag{ 4 }$$

where ρ_0,d is the midplane density, Φ_0,d is the midplane potential, and c_s,d is the isothermal sound speed in the disk, ${c}_{{\rm{s}},{\rm{d}}}=\sqrt{{k}_{{\rm{B}}}{T}_{\mathrm{disk}}/\mu {m}_{{\rm{p}}}}$. When converting between mass density ρ and number density n, we take μ = 0.6 throughout the calculation, as appropriate for ionized gas with solar metallicity. The midplane density (at a given radius) is calculated by requiring the integrated density profile to equal the surface density, so

$$\begin{eqnarray}&&{\rho }_{{\rm{0,d}}}=\displaystyle \frac{{\rm{\Sigma }}}{{\int }_{-\infty }^{\infty }\exp \left(-\tfrac{{\rm{\Phi }}-{{\rm{\Phi }}}_{{\rm{0,d}}}}{{c}_{{\rm{s,d}}}^{2}}\right)}.\end{eqnarray} \tag{ 5 }$$

The disk gas is also in radial equilibrium with the static potential. Velocities are set by first calculating the tangential acceleration at a given radius due to the gravitational potential, then correcting this acceleration for the radial pressure gradient, i.e.,

$$\begin{eqnarray}&&{a}_{\phi }(r,z)=-{\rm{\nabla }}{\rm{\Phi }}(r,z)+\displaystyle \frac{1}{\rho }\displaystyle \frac{{dP}}{{dr}},\end{eqnarray} \tag{ 6 }$$

where $P=\rho {c}_{{\rm{s}}}^{2}/\gamma$ is the gas pressure. We assume an adiabatic index γ = 5/3 throughout the simulations. The disk gas is artificially truncated via an exponential ramp down in surface density beyond a radius of R = 4.5 kpc, in order to reduce possible boundary effects from the simulation box.

2.1.3. Halo Gas

The isothermal gas disk is embedded in an adiabatic hydrostatic halo. To set the initial conditions for the halo gas, we first determine the density profile as a function of spherical radius, which for an adiabatic gas in hydrostatic equilibrium with a spherical potential is given by

$$\begin{eqnarray}&&\rho (r)={\rho }_{0,{\rm{h}}}{\left[1+(\gamma -1)\displaystyle \frac{{\rm{\Phi }}-{{\rm{\Phi }}}_{0,{\rm{h}}}}{{c}_{{\rm{s}},{\rm{h}}}^{2}}\right]}^{\displaystyle \frac{1}{\gamma -1}}.\end{eqnarray} \tag{ 7 }$$

In this case, we normalize the profile by manually setting the density, ρ_0,h, at a radius of 100 kpc. We assume ${\rho }_{0,{\rm{h}}}=3\times {10}^{3}\,{M}_{\odot }\,{\mathrm{kpc}}^{-3}$, or a number density of around $n\approx {10}^{-3.5}\,{\mathrm{cm}}^{-3}$. Similarly, Φ_0,h is the gravitational potential at the cooling radius. Unlike the disk, in setting the halo gas distribution the potential is taken to be spherically symmetric (in other words, we evaluate the disk component of Φ at r = 0 for all z). Finally, the sound speed in Equation (7) is set by assuming a temperature of T = 10⁶ K at 100 kpc. The normalization of the number density and temperature we have chosen for the halo gas may be more appropriate for a slightly larger galaxy. However, these values for the M82 halo at 100 kpc are highly uncertain given its location within the larger M81 group. In practice, the halo gas gets blown out of the domain at early times in our simulations, so the exact normalization does not affect our results. Having set the density, we then set the pressure of the gas adiabatically using P = Kρ^γ, where $K={c}_{{\rm{s}},{\rm{h}}}^{2}{\rho }_{0,{\rm{h}}}^{1-\gamma }/\gamma$.

The density and internal energy from the halo profile are added to every cell in the simulation volume, including those in the disk, ensuring a smooth transition between disk and halo. Figure 4 shows the resulting one-dimensional density and temperature profiles as a function of radius along the z-axis and in the disk midplane. We have determined that these initial conditions are stable for over 1 Gyr, well beyond the timescale of our galactic wind simulations. Shear forces at the interface between the rotating disk gas and the static gaseous halo do generate slight turbulence, but despite these perturbations the disk remains globally stable over long periods.

Zoom In Zoom Out Reset image size

2.2. Feedback Model

In this paper we will focus on the first, simplest feedback model that is implemented in the CGOLS suite. The series of simulations is designed to be systematic, building complexity in the feedback model from one simulation to the next, in order to better constrain our analytic models of feedback and winds. The first wind model we test is therefore that of a nuclear starburst driven by continuous star formation. The original analytic model for such a wind was developed by Chevalier & Clegg (1985, hereafter CC85) in order to explain the wind observed in the nearby starburst galaxy M82. The CC85 model assumes adiabatic hydrodynamics and neglects the effects of gravity. We note that there have been several theoretical improvements to the CC85 model that include additional physics (Wang 1995; Silich et al. 2004; Bustard et al. 2016; Thompson et al. 2016), which are addressed in our companion paper.

In the CC85 model, a constant source of mass and energy are deposited in a spherical volume with radius R. The mass input rate $\dot{M}$ can be quantified as a function of the star formation rate, ${\dot{M}}_{\mathrm{SFR}}$,

$$\begin{eqnarray}&&\dot{M}=\beta {\dot{M}}_{\mathrm{SFR}},\end{eqnarray} \tag{ 8 }$$

where β is known as the mass-loading factor of the wind. Similarly, the energy input rate $\dot{E}$ is quantified as a fraction of the energy input by supernovae ${\dot{E}}_{\mathrm{SN}}$,

$$\begin{eqnarray}&&\dot{E}=\alpha {\dot{E}}_{\mathrm{SN}},\end{eqnarray} \tag{ 9 }$$

where the energy-loading factor, α, accounts for radiative losses in the ISM. With the assumption that each supernova releases 10⁵¹ erg of energy, that there is one supernova per $100\,{M}_{\odot }$ of star formation, and that the star formation rate is in units of M_⊙ yr⁻¹, the energy loading can be rewritten as a function of the star formation rate,

$$\begin{eqnarray}&&\dot{E}=3\times {10}^{41}\,\mathrm{erg}\,{{\rm{s}}}^{-1}\,\alpha {\dot{M}}_{\mathrm{SFR}}.\end{eqnarray} \tag{ 10 }$$

To replicate this model in the initial CGOLS simulations, we continuously inject mass and energy into a spherical volume in the center of the galaxy disk. We use a radius of R = 300 pc for the gain region. Cells along the edge of the sphere are weighted so as to recover the correct total injection rate at low resolution and alleviate grid effects (we find the weighting unnecessary in our higher-resolution simulations). Over the course of the simulation, we vary α, β, and the star formation rate (SFR) in order to probe different wind regimes. Although radiative losses will not be discussed in this work, different mass- and energy-loading factors lead to different theoretical predictions regarding the importance of radiative cooling in the wind. Varying α, β, and SFR allows us to test different theoretical models within a single simulation. We refer to the feedback model described here as "central" in Table 1. Other feedback models are discussed in additional CGOLS papers, including a clustered model that distributes the mass and energy nonisotropically throughout the central disk region (Schneider et al. 2018).

We start the simulations by running for 5 Myr with no feedback, primarily to allow cells at the disk–halo interface to equilibrate in radiative simulations. At 5 Myr, we begin the mass and energy injection with $\dot{M}=1.5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and $\dot{E}=1.5\times {10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, corresponding to β = 0.3 and α = 1.0 with an assumed star formation rate of ${\dot{M}}_{\mathrm{SFR}}=5\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. We refer to this throughout the paper as the low mass-loading state. Immediately after starting the feedback, we ramp up to a higher mass-loading state, with $\dot{M}=12\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$ and $\dot{E}=5.4\times {10}^{42}\,\mathrm{erg}\,{{\rm{s}}}^{-1}$, corresponding to β = 0.6 and α = 0.9 with an assumed star formation rate of ${\dot{M}}_{\mathrm{SFR}}=20\,{M}_{\odot }\,{\mathrm{yr}}^{-1}$. The ramping is a linear function of time, e.g.,

$$\begin{eqnarray}&&\dot{M}={M}_{1}+\displaystyle \frac{t}{{t}_{\mathrm{ramp}}}({M}_{2}-{M}_{1})\end{eqnarray} \tag{ 11 }$$

for the mass loading, where M₁ is the low mass-loading state, M₂ is the high mass-loading state, t_ramp = 5 Myr, and t is the total time since the feedback began. We then keep the simulation in the high mass-loading state for 30 Myr, before ramping back down to the low mass-loading state, e.g.,

$$\begin{eqnarray}&&\dot{M}={M}_{1}+\displaystyle \frac{t-40\,\mathrm{Myr}}{{t}_{\mathrm{ramp}}}({M}_{1}-{M}_{2}),\end{eqnarray} \tag{ 12 }$$

again over the course of 5 Myr. The remaining 30 Myr is run in the low mass-loading state, for a total simulation run time of 75 Myr. The mass and energy loading as a function of time are plotted in Figure 5. In addition to sampling an interesting region of parameter space, these α, β, and SFR values are also physically motivated by the starburst itself. At the simulation start, we assume that some outside event (in this case, an interaction with nearby M81) has driven a large amount of gas to the center of the galaxy, triggering a burst of star formation. The mass loading may be higher at early stages in the burst because the supernova remnants are interacting with ISM gas that can add additional mass to the hot wind. Correspondingly, α is lower at early times as a result of radiative losses associated with these interactions (Kim et al. 2017). As the burst matures, the star formation rate decreases, β drops, and α increases as much of the ISM has been cleared out. The late-time low mass-loading parameters are directly motivated by constraints from X-ray observations of the hot wind in M82 today (Strickland & Heckman 2009).

Zoom In Zoom Out Reset image size

2.3. Hydrodynamical Model

In this section, we present scientific results from the first simulation carried out in the CGOLS suite, the adiabatic simulation A-2048. The most straightforward of the models, it consists of a centrally driven wind and neglects the effects of radiative cooling. Thus, this model serves as a test of the original CC85 analytic outflow model within the context of a disk and halo and provides a basis of comparison for the more realistic models that are presented in other papers in this series.

In Figure 6, we show density and temperature slices through the x–z plane of the simulation at a series of characteristic times: t = 10 Myr, when the initial superbubble is moving out into the CGM; t = 25 Myr, about 15 Myr after the high mass-loading wind state has commenced; t = 50 Myr, when the ramp down has occurred and the low mass-loading outflow solution is propagating outward; and t = 60 Myr, about 15 Myr after the simulation reverts to the low mass-loading wind state. Relevant features of the outflow at each snapshot time include:

Zoom In Zoom Out Reset image size

• 1.  
  t = 10 Myr: The superbubble created by the mass and energy input at the galaxy's center is propagating out into the CGM. The forward shock, contact discontinuity, and reverse shock are visible in both the density and temperature slice. Behind the reverse shock, the flow is smooth, and the beginnings of a pressure-confined conical outflow region can be seen. The isothermal disk remains relatively undisturbed at this time, though its central regions were blown out by the mass and energy injection immediately following the start of feedback. The outflow is just beginning to interact with the surface of the disk.
• 2.  
  t = 25 Myr: At this point, the high mass-loading outflow has stabilized in a biconical region centered on the z-axis. The free-wind cone is surrounded by a region of higher temperature and density caused by interactions between the outflow and disk. The wind from the central outflow has begun to eat away at the high-density gas near the center of the disk, driving turbulence along the interface, and the shock from the initial bubble has propagated through the disk, raising the temperature as compared to the initial conditions. This configuration remains more or less constant for the following 15 Myr of evolution in the high mass-loading state.
• 3.  
  t = 50 Myr: The simulation is now 5 Myr into the lower mass-loading state. While the outflow solution remained smooth during the ramp down from 40–45 Myr, the period t = 45–50 Myr is characterized by a second outward-moving shock caused by the new outflow solution overtaking the old. The lower mass-loading state and higher α/β ratio results in a significantly higher terminal velocity, as expected from the analytic estimate,

  $$\begin{eqnarray}&&{v}_{\mathrm{term}}=\sqrt{2\dot{E}/\dot{M}},\end{eqnarray} \tag{ 13 }$$

  which gives v_term = 1776 km s⁻¹ and v_term = 1191 km s⁻¹ for the low and high mass-loading states, respectively. The turbulence at the interface has now increased considerably, and the resulting high pressure region combined with the lower pressure of the new outflow state is forcing the undisturbed outflow into an even narrower cone.
• 4.  
  t = 60 Myr: After 15 Myr in the low mass-loading state, the turbulent high pressure region from the interaction between the outflow and disk has overtaken most of the simulation volume. The undisturbed outflow now resembles a cylinder centered along the z-axis, and the disk itself has puffed up relative to the initial conditions as a result of the increased temperature and shocks that have now crossed it multiple times. Large clouds of higher-density disk gas continue to be ablated from the disk–wind interface and are quickly shredded and carried away in the flow. The bow shocks from these interactions result in significantly higher temperature regions in the interface region than in the undisturbed outflow cone.

3.1. Comparison with the Analytic Model

Within the region of relatively undisturbed flow, we find an excellent fit between our simulations and the CC85 model. The CC85 model gives a solution for the spherically expanding flow in terms of the Mach number, which transitions from subsonic to supersonic at the edge of the mass and energy input region, R. The solutions for other physical variables of interest, such as density, velocity, and pressure, can be found as a function of radius by numerically integrating the spherical hydrodynamic equations with mass and energy source terms, e.g.,

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{d}{{dr}}(\rho {v}_{r}{r}^{2})=\displaystyle \frac{\dot{M}}{V},\end{eqnarray} \tag{ 14 }$$

$$\begin{eqnarray}&&\rho {v}_{r}\displaystyle \frac{{{dv}}_{r}}{{dr}}=-\displaystyle \frac{{dP}}{{dr}}-\displaystyle \frac{\dot{M}}{V}{v}_{r}\end{eqnarray} \tag{ 15 }$$

$$\begin{eqnarray}&&\displaystyle \frac{1}{{r}^{2}}\displaystyle \frac{d}{{dr}}\left[\rho {v}_{r}{r}^{2}\left(\displaystyle \frac{1}{2}{v}_{r}^{2}+\displaystyle \frac{\gamma }{\gamma -1}\displaystyle \frac{P}{\rho }\right)\right]=\displaystyle \frac{\dot{E}}{V},\end{eqnarray} \tag{ 16 }$$

where v_r is the radial velocity, $V=\tfrac{4}{3}\pi {R}^{3}$ is the volume within the injection region, and all other symbols are as previously defined.

In Figure 7, we compare the analytic solution predicted by the CC85 model to the results from the simulation at t = 25 Myr (in the high mass-loading state) and in Figure 8 at 60 Myr (in the low mass-loading state). Panels show the density, radial velocity, and pressure, with the analytic solution plotted as a black line. In order to focus the results on the outflow, only simulation data within a biconical region with an opening angle of ΔΩ = 60° centered on the z-axis are included in this analysis. Simulation data are binned as a function of spherical radius into 93 evenly spaced bins with Δr = 0.125 kpc.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

At early times, the simulation data within the outflow cone match the analytic solution nearly exactly. The maximum velocity is slower by about 40 km s⁻¹, which can be attributed to the effect of the gravitational potential that is not taken into account in the analytic model. At later times, the results from the simulation appear more complex. The outflow cone during the low mass-loading state is narrower and less time-stable, effectively "wobbling" around as the pressure balance with the turbulent interface region changes. As a result, much of the volume within the measured biconical region is now filled with turbulent interface material. The turbulent region is characterized by gas that has higher density and lower velocity, leading to an overall increase in pressure and temperature, as seen at larger radii in Figure 8. The shocks created by slowing the hot wind gas as it interacts with clouds of denser, cooler gas can explain this increased temperature. For example, at r = 10 kpc, the mean radial velocity of the interface material measured within the ΔΩ = 60° bicone is v_r = 1610 km s⁻¹, approximately 150 km s⁻¹ less than the terminal velocity predicted by the analytic model. Slowing the gas by this amount would increase its temperature by ΔT ∼ 9 × 10⁵ K, and indeed, we measure a mean temperature for the gas at r = 10 kpc of T = 1.05 × 10⁶ K, while the analytic expectation for the outflow temperature at r = 10 kpc is T ∼ 10⁵ K. This slowing of the outflowing material and increased density and temperature should also affect the observable signature of soft X-ray emission, which we investigate in the following section.

Many studies of the soft X-ray emission have been conducted for M82 (e.g., Griffiths et al. 2000; Strickland et al. 2004; Li & Wang 2013), one of the best-resolved starburst galaxies in the local universe and a unique target for investigating the spatial distribution of diffuse soft X-ray emission. We compare the soft X-ray emission from our adiabatic simulation to observations of the emission in M82, in order to both validate our model and provide additional insight as to the source of the observed emission. We focus our comparison on a snapshot from the simulation at t = 60 Myr, as the later stages of the simulation represent the most relevant comparison to the present-day M82 system. The mass and energy loading during the late-time evolution of the simulation are set to match those of M82 by construction, based on estimates of α and β made by studying the diffuse hard X-ray emission in the nuclear region (Strickland & Heckman 2009). However, those constraints do not necessarily ensure a good match to the soft X-rays.

In Figure 9, we show a map of the soft X-ray surface brightness at t = 60 Myr. The map was made by integrating the estimated soft X-ray emission along the y-axis, where only cells with temperatures between 0.3 and 2.0 keV were included in the integral, i.e.,

$$\begin{eqnarray}&&{S}_{{\rm{X}}}=\displaystyle \frac{1}{4\pi }{\int }_{0}^{{L}_{y}}{n}^{2}{\rm{\Lambda }}(T){dy}\ [\mathrm{erg}\,{{\rm{s}}}^{-1}\,{\mathrm{cm}}^{-2}\,{\mathrm{str}}^{-1}].\end{eqnarray} \tag{ 17 }$$

We convert from temperature in a given cell to energy by assuming $E=\tfrac{3}{2}{k}_{{\rm{B}}}T$, as appropriate for a fully ionized plasma. The cooling curve used to estimate the emission, Λ(T), is an analytic fit to a solar-metallicity collisional ionization equilibrium curve calculated with Cloudy (Ferland et al. 2013) and is described in Appendix A.3. The line-of-sight integral gives a brightness in erg s⁻¹ cm⁻², which we convert to surface brightness by dividing by 4π steradian, yielding the values shown in the map in Figure 9. The 0.3 − 2.0 keV band was chosen for comparison to the Chandra X-ray telescope ACIS data.

Zoom In Zoom Out Reset image size

Overall, the simulated X-ray data compare favorably to the observations in both spatial extent and luminosity (see, e.g., the surface brightness map in Figure 2(h) of Strickland et al. 2004). The integrated soft X-ray luminosity across the entire simulation volume is L_X = 1.9 × 10⁴⁰ erg s⁻¹. Strickland et al. (2004) calculate a total soft X-ray luminosity of L_X = 4.3 × 10⁴⁰ erg s⁻¹, of which L_X = 2.3 × 10⁴⁰ erg s⁻¹ comes from the nuclear region, a circular aperture 1 kpc in radius. Given that this simulation is adiabatic, the extent to which the soft X-ray emission matches the observations should perhaps be surprising, but we do note that in this region of α and β parameter space, we do not expect significant radiative losses from the outflow.

For a more quantitative look at how the spatial extent of the soft X-rays compares with the data, we plot in Figure 10 the vertical surface brightness profile of the map shown in Figure 9. We integrate the surface brightness as calculated in Equation (17) in slices with Δx = 10 kpc and Δz = 0.15625, and then divide over the area to keep the units of erg s⁻¹ cm⁻² arcmin⁻² (with a conversion from steradians to arcminutes). We have also plotted the best-fit exponential profile from Strickland et al. (2004), which was calculated in the same manner. To facilitate comparison with the observations, we have converted Chandra counts in the 0.3 − 2.0 keV band to energies assuming 1 photon ∼ 2 × 10⁻⁹ erg (see, e.g., Section 4.3.2 of Strickland et al. 2004).

Zoom In Zoom Out Reset image size

Again, we find a good match between the simulated data and the observations. The extended filamentary structures seen in the simulation lead to a surface brightness profile that is less peaked than what would be produced by a pure CC85 outflow. Figure 9 shows that the soft X-rays in the simulation are highly collimated by the disk. This degree of collimation produced in the simulation is likely not realistic, given that the adiabatic nature of the disk leads it to be much puffier at late times than would be expected if it were allowed to cool radiatively.

The high-resolution simulations in the CGOLS suite are orders of magnitude larger than those that have been used in other contexts to study galactic winds (see Figure 1). In this section, we investigate why this level of resolution is necessary. In order to test for convergence of measured quantities, all of the Δx ≈ 5 pc simulations in the CGOLS suite are additionally carried out at both Δx ≈ 10 pc and Δx ≈ 20 pc resolution. We note that these comparison runs are still large simulations—a Δx ≈ 20 pc resolution simulation has a total of 512 × 512 × 1024 ∼ 250 million cells. Nevertheless, we find that the overall nature of the simulations changes significantly as a result of the 2× improvement in resolution between the A-512 and A-1024 models.

At early times, the results from all three simulations are converged. The radial outflow data presented in Figure 7 for the outflow cone are identical between all three simulations. Other features of the simulations at early times are converged as well, for example, the shock speed of the superbubble as it travels through the ambient CGM and leaves the simulation volume. At later times, we find that the large-scale features remain converged—there is little difference in the radial data plots shown in Figure 8, for example. However, the degree of turbulence generated in the interface region is quite different for the low-resolution A-512 simulation compared to higher resolutions. As an example, we show in Figure 11 the 60 Myr temperature slice for the A-512 simulation, which can be directly compared with the slice shown for the fiducial A-2048 simulation in Figure 6.

Zoom In Zoom Out Reset image size

We can better quantify the effects of resolution by estimating the degree to which turbulent gas plays a role in each simulation. To do this, we calculate the rms velocity as a function of z-height above the disk midplane. In the nonturbulent case, the velocity at any given point in the outflow is expected to be radial, so to determine the contribution of turbulent gas, we calculate for every cell the magnitude difference between the volume-weighted total velocity, $v=\sqrt{{v}_{x}^{2}+{v}_{y}^{2}+{v}_{z}^{2}}$, and the volume-weighted radial velocity, v_r = (v_x x + v_y y + v_z z)/r. We then determine the rms of this velocity difference,

$$\begin{eqnarray}&&{\rm{\Delta }}{v}_{\mathrm{rms}}=\sqrt{\langle {(v-{v}_{r})}^{2}\rangle },\end{eqnarray} \tag{ 18 }$$

for all cells in the simulation as a function of z-height above the disk. The resulting profiles are plotted in Figure 12. At low z, Δv_rms primarily shows the average rotation velocity of gas in the disk. Starting at z ∼ 1 kpc, however, the increase in turbulent velocities at higher resolution is clear and is particularly notable for the highest-resolution simulation at large z. A mass-weighted version of Figure 12 shows a similar trend, indicating that both the dense gas and the volume-filling gas are turbulent primarily at high resolution.

Zoom In Zoom Out Reset image size

Because the outflow velocities are expected to primarily be radial, we interpret Δv_rms as a measure of the turbulence in the outflow. (We have checked that the values of Δv_rms cannot be explained by residual rotational motions of the disk gas.) In these simulations, this turbulence is in part generated by shear instabilities between the radially moving hot wind and the rotating disk gas. Development of the Kelvin–Helmholtz instability at unresolved interfaces is known to be affected by resolution (e.g., Robertson et al. 2010), and the growth rate of the instability is fastest for the smallest wavelength perturbations. As a result, perturbations at the disk–wind interface in the higher-resolution simulations will more quickly generate instability on smaller scales. These small eddies are effective in lifting clouds of gas out of the disk and into the outflow, where they contribute to the measured Δv_rms. Given the relatively short crossing time for the hot gas to traverse the disk (∼3–5 Myr), a resolution of 20 pc appears insufficient to allow development of the instability at a level that significantly affects turbulent velocities in the outflow.

In addition to the turbulent velocities, we might also expect the clumpiness of the gas in the outflow to be affected by the resolution of the simulation. In order to check this, we additionally plot in Figure 13 the clumping factor of the gas as a function of z-height for the three different resolutions. We bin all of the data in the simulations in slices of Δz = 0.1 kpc and calculate the dimensionless clumping factor for each slice, $\langle {n}^{2}\rangle /\langle n{\rangle }^{2}$. Perhaps surprisingly, we see relatively little difference in the clumpiness between the three simulations. As the density slices in Figure 6 lead us to expect, all three simulations show the most clumping at heights close to the disk, and progressively less as z increases. There may be a slight trend toward increased variability as a function of resolution, and the two higher-resolution simulations show significantly more clumpiness at high latitudes than model A-512.

Zoom In Zoom Out Reset image size

The turbulence in the outflow is supersonic and produces shocks. Thus, the increase in turbulent structure as a function of simulation resolution has observable consequences for the spatial distribution of the soft X-ray emitting gas. In Figure 14, we show the soft X-ray surface brightness maps for the A-512, A-1024, and A-2048 simulations. Clearly, the overall nature of the outflow has not yet converged at Δx = 20 pc resolution (the leftmost map). We note that this is despite the fact that the total size of the simulation volume is much larger than the resolution of a single cell, and the region where the mass and energy is being deposited is also resolved with many cells, given that R = 300 pc. Large-scale shocks propagating throughout the domain tend to dominate the features seen in the X-ray emission map and are also the dominant feature in temperature and density projections (not shown). At the higher resolutions, both the density and temperature projections show a similar filamentary structure, though not as obviously as in the soft X-ray plot. The filaments arise as a result of the shredding and ablation of the denser gas clouds that have been lofted into the outflow from the disk (in movies⁵ , we can watch these clouds as they move outward and more easily observe the expansion and shredding). Many of these clouds come from the inner region of the disk near the edge of the injection zone, but some additional dense gas is contributed from the shear instabilities at the disk–outflow interface, particularly at later times.

Zoom In Zoom Out Reset image size

Whether the total soft X-ray luminosity is converged even in the highest-resolution simulation (Model A-2048) is not clear. For the two higher-resolution simulations, L_X tends to increase with resolution. In simulation A-512, the nuclear region and large-scale shock features contribute significantly to the soft X-ray luminosity, such that we find the largest total luminosity in the simulation with the lowest resolution, but the relationship is not monotonic: L_X = 2.9 × 10⁴⁰ erg s⁻¹, L_X = 1.4 × 10⁴⁰ erg s⁻¹, and L_X = 1.9 × 10⁴⁰ erg s⁻¹ for A-512, A-1024, and A-2048, respectively. However, despite the total luminosity in simulation A-512 being the closest to the observed luminosity, the z-profile of the emission in the low-resolution simulation is much more steeply peaked than is observed. We do not wish to overinterpret these differences in total luminosity, given the uncertainties in the estimated emission caused by our analytic fit to the cooling curve and the simplifying assumption that only cells with volume-average temperature between 0.3 and 2.0 keV contribute to the X-ray emission. At the very least, we are confident that the match to the morphology of the observed soft X-ray brightness profile is improved at higher resolution.

In this work, we have introduced the CGOLS suite, a set of high-resolution isolated-galaxy simulations with feedback designed to generate a galactic wind. The volume of the simulation boxes (10 × 10 × 20 kpc) allows us to study the evolution of the winds as they escape the galaxy and begin to interact with the CGM, while the resolution (Δx < 5 pc across the entire volume) is sufficient to capture the small-scale hydrodynamics and thermal instabilities that contribute to mixing and potential formation of multiphase gas.

In addition to a detailed presentation of the features common to all simulations in the suite, we have focused the results of this paper on the first of the CGOLS simulations, an adiabatic simulation of a galaxy with a nuclear starburst modeled after M82 (model A-2048). At early times, the analytic model of a nuclear starburst described by Chevalier & Clegg (1985) fits the data generated by the simulation well, within a biconical region with an opening angle of ΔΩ = 60°. Outside this region and at late times, significant turbulence is generated in the interface region between the disk gas and the central hot outflow. This region is dominated by gas with densities and temperatures up to an order of magnitude higher than those predicted by the analytic model.

Even with the simplified adiabatic physics of the A-2048 simulation, we calculate a total luminosity in soft X-ray emission that is only a factor of 2 lower than the observed soft X-ray emission in the starburst galaxy M82 (see Figure 10). The turbulent interface region is a key feature in correctly generating the shape of the X-ray profile as a function of height above the disk, given that interactions between the hot, fast wind and slower, denser clouds of gas generate shocks and mixing at large radii that contribute significantly to the X-ray emission.

This paper is the first in a series that will systematically investigate the nature of galactic winds through large-scale global galaxy simulations. Results from additional simulations in the CGOLS suite that include radiative cooling and probe the multiphase nature of the outflow are presented in a companion paper (Schneider et al. 2018).

We acknowledge inspiration from the Simons Symposia on Galactic Superwinds: Beyond Phenomenology. This research used resources of the Oak Ridge Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC05-00OR22725, via an award of computing time on the Titan system provided by the INCITE program (AST125). E.E.S. is supported by NASA through Hubble Fellowship grant No. HF-51397.001-A awarded by the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., for NASA, under contract NAS 5-26555. B.E.R. is partially supported by NASA grant 80NSSC18K0563, contract NNG16PJ25C, and program HST-GO-14747.

Software: Cholla (Schneider & Robertson 2015), matplotlib (Hunter 2007), numpy (Van Der Walt et al. 2011), hdf5 (The HDF Group 1997-NNNN), Cloudy (Ferland et al. 2013).

A.1. Hydrodynamics Solver

This work used a hydrodynamics integration algorithm not previously incorporated in Cholla, so we document it here for completeness. The solver is based on the description provided in Stone & Gardiner (2009) of the "Van Leer" integration scheme for the Athena code, a predictor–corrector scheme modeled in turn after the MUSCL-Hancock scheme (see, e.g., Toro 2009). This second-order integrator updates the vector of conserved variables

$$\begin{eqnarray}&&{\boldsymbol{u}}={[\rho ,\rho u,\rho v,\rho w,E]}^{T}\end{eqnarray} \tag{ 19 }$$

from time t = n to time t = n + 1 using fluxes estimated at the half time step $t=n+\tfrac{1}{2}$. Here, ρ is the mass density; u, v, and w are the three components of the velocity vector; and E is the total energy density. In three dimensions, the update for a cell (i, j, k), can be written

$$\begin{eqnarray}&&{{\boldsymbol{u}}}_{(i,j,k)}^{n+1}={{\boldsymbol{u}}}^{n}+\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}x}[{{\boldsymbol{F}}}_{\left(i,-,\tfrac{1}{2},,,j,,,k\right)}^{n+\tfrac{1}{2}}-{{\boldsymbol{F}}}_{\left(i,+,\tfrac{1}{2},,,j,,,k\right)}^{n+\tfrac{1}{2}}]\\ &&+\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}y}[{{\boldsymbol{G}}}_{\left(i,,,j,-,\tfrac{1}{2},k\right)}^{n+\tfrac{1}{2}}-{{\boldsymbol{G}}}_{\left(i,,,j,+,\tfrac{1}{2},,,k\right)}^{n+\tfrac{1}{2}}]+\displaystyle \frac{{\rm{\Delta }}t}{{\rm{\Delta }}z}[{{\boldsymbol{H}}}_{\left(i,,,j,,,k,-,\tfrac{1}{2}\right)}^{n+\tfrac{1}{2}}-{{\boldsymbol{G}}}_{\left(i,,,j,,,k,+,\tfrac{1}{2}\right)}^{n+\tfrac{1}{2}}].\end{eqnarray} \tag{ 20 }$$

Here, ${{\boldsymbol{F}}}^{n+\tfrac{1}{2}}$, ${{\boldsymbol{G}}}^{n+\tfrac{1}{2}}$, and ${{\boldsymbol{H}}}^{n+\tfrac{1}{2}}$ are estimates of the flux between cells at the half time step; Δt is the hydrodynamic time step; and Δx, Δy, and Δz are the cell sizes in each dimension. The basic steps of the three-dimensional algorithm are as follows:

• 1.  
  Construct first-order accurate fluxes at each cell interface using the cell-averaged values of the conserved variables and the Riemann solver of choice (in this work, we use the HLLC solver). To calculate the first-order flux ${{\boldsymbol{F}}}_{i+\tfrac{1}{2}}^{* }$ at the interface between cells i and i + 1, solve the Riemann problem beginning with left state values ${{\boldsymbol{W}}}_{L}=({\rho }_{i}$, v_i, p_i) and right state values ${{\boldsymbol{W}}}_{R}=({\rho }_{i}+1$, v_i+1, p_i+1).
• 2.  
  Use these first-order accurate fluxes to update the conserved variables in each cell by a half time step, $\tfrac{{\rm{\Delta }}t}{2}$. Apply the update from Equation (20) with ${\rm{\Delta }}t=\tfrac{{\rm{\Delta }}t}{2}$ and using ${\boldsymbol{F}}* ,{\boldsymbol{G}}*$, and ${\boldsymbol{H}}*$.
• 3.  
  Use the half-step values of the conserved variables, ${{\boldsymbol{u}}}^{n+\tfrac{1}{2}}$, to construct time-centered estimates of the conserved variables at the cell edges, ${{\boldsymbol{U}}}_{L}^{n+\tfrac{1}{2}}$ and ${{\boldsymbol{U}}}_{R}^{n+\tfrac{1}{2}}$, using the interpolation method of choice (in this work, we use a piecewise linear reconstruction with limiting carried out in the characteristic variables). Note that unlike PPM (Colella & Woodward 1984), this algorithm does not require a time evolution of the reconstructed interface values.
• 4.  
  Solve a second set of Riemann problems using the time-centered spatially reconstructed variables, ${{\boldsymbol{U}}}_{L}^{n+\tfrac{1}{2}}$ and ${{\boldsymbol{U}}}_{R}^{n+\tfrac{1}{2}}$, to calculate the time-centered second-order fluxes, ${{\boldsymbol{F}}}^{n+\tfrac{1}{2}}$, etc.
• 5.  
  Use Equation (20) with the full time step and the second-order fluxes to calculate the new values of the conserved variables at time n + 1.

In addition to being relatively simple and robust, this algorithm provides a set of first-order fluxes, ${\boldsymbol{F}}* ,\,{\boldsymbol{G}}\,*$, and ${\boldsymbol{H}}*$. When used in Equation (20), these fluxes have the advantage of being positive definite—that is, in the pure hydrodynamic case, they will not produce negative densities or energies after the cell update. Given the rather extreme nature of the calculations presented in this work, we have found that the full integration scheme can still sometimes produce negative densities or total energies for a given cell when used with higher-order reconstruction methods, particularly in simulations where the gas is also allowed to cool radiatively. In that case, we use the first-order fluxes to perform the complete update, only for the affected cell. In practice, this tends to affect of the order of 1 × 10⁻¹⁰ cells in the calculation on any given update and was never necessary in the adiabatic simulation presented in this paper.

A.2. Static Gravity Implementation

In order to carry out the simulations presented in this work, we added a static gravity module to the Cholla code. The time-averaged gravitational source terms are coupled to the hydrodynamics module using an operator-split cell update, meaning that first each cell is updated using the hydrodynamic fluxes according to Equation (20), then the gravitational source terms are applied to the new momentum and energy. The source terms for momentum and energy are defined as

$$\begin{eqnarray}{{\boldsymbol{S}}}_{m} & = & \rho {\rm{\nabla }}{\rm{\Phi }}=\displaystyle \frac{{\rm{\Delta }}t}{2}{\boldsymbol{g}}({\rho }^{n}+{\rho }^{n+1}),\quad {{\boldsymbol{S}}}_{E}=\rho {\boldsymbol{v}}\cdot {\rm{\nabla }}{\rm{\Phi }}\\ & = & \displaystyle \frac{{\rm{\Delta }}t}{4}{\boldsymbol{g}}\cdot ({{\boldsymbol{v}}}^{n}+{{\boldsymbol{v}}}^{n+1})({\rho }^{n}+{\rho }^{n+1}),\end{eqnarray} \tag{ 21 }$$

where ${\boldsymbol{g}}$ is the gravitational acceleration calculated from the potential Φ, ρⁿ and ρⁿ⁺¹ are the density before and after the hydrodynamic update, and ${{\boldsymbol{v}}}^{n}$ and ${{\boldsymbol{v}}}^{n+1}$ are the velocity vectors before and after the hydrodynamic update.

A.3. Cooling Implementation

Although we do not use radiative cooling in the simulation described in this paper, we use it for the rest of the simulations in the CGOLS suite, so we include our implementation here for completeness. We use an operator-split approach to account for radiative losses in our simulations, as described in Appendix C of Schneider & Robertson (2017). The only difference is in the calculation of the cooling function, Λ. Rather than using Cloudy tables as described in that work, here we use a simple analytic function that is a piecewise parabolic fit to a solar-metallicity collisional ionization equilibrium cooling curve computed with Cloudy. The functional form of the cooling curve is

$$\begin{eqnarray}{\rm{\Lambda }} & = & 0,\quad {10}^{4}\lt T\\ {\rm{\Lambda }} & = & {10}^{-1.3{(T-5.25)}^{2}-21.25},\quad {10}^{4}\leqslant T\lt {10}^{5.9}\\ {\rm{\Lambda }} & = & {10}^{0.7{(T-7.1)}^{2}-22.8},\quad {10}^{5.9}\leqslant T\lt {10}^{7.4}\\ {\rm{\Lambda }} & = & {10}^{0.45T-26.065},\quad {10}^{7.4}\leqslant T,\end{eqnarray} \tag{ 22 }$$

where T is in logarithmic units. The cooling curve and fit are shown in Figure 15. The simulations are run with a temperature floor of T = 10⁴ K, so we do not take into account cooling below this threshold, nor do we account for any radiative heating processes. This cooling curve is also used to calculate the soft X-ray emission maps presented in Section 4.

Zoom In Zoom Out Reset image size