Galactic winds - How to launch galactic outflows in typical Lyman-break galaxies
We perform hydrodynamical simulations of a young galactic disc embedded in a hot gaseous halo using parameters typical for Lyman break galaxies (LBGs). We take into account the (static) gravitational potentials due to a dark matter halo, a stellar bulge and a disc of stars and gas. Star formation is treated by a local Kennicutt-Schmidt law. We simplify the structure of the interstellar medium by restricting the computational domain to a 25th of the full azimuthal angle, effectively assuming large-scale axisymmetry and neglecting any effects of spiral structure, and focus on the large-scale ISM drivers, the superbubbles. Supernovae are triggered randomly and have preset event sizes of several tens to hundreds. We further investigate different halo gas pressures and energy injection methods. Many of our simulated galaxies, but not all, develop bipolar outflows. We characterise the strength of the outflow by mass and energy outflow rates, and investigate the effect of changes to the details of the model. We find that supernovae are more effective if comprised into larger superbubbles. The weight and the pressure of the halo gas is able to quench galactic outflows. The wind emerges from a series of superbubbles in regions where a critical star formation density is exceeded. The superbubbles expand into the gaseous halo at slightly supersonic speed, producing radiative shock waves with similar characteristics as the absorptions systems observed around LBGs.
keywords:methods: numerical – galaxies: high-redshift – galaxies: evolution – galaxies: haloes – galaxies: ISM – ISM: bubbles
Galactic winds are commonly characterised by supersonic, biconically shaped outflows perpendicular to the midplane of a galactic disc, energetic enough to carry significant amounts of gas far away from their host system, and in some cases even into the intergalactic medium (IGM, for a review see Veilleux, Cecil & Bland-Hawthorn (2005)). Observations have revealed evidence for supersonic outflows in the spectra of nearby ultra-luminous infrared galaxies (ULIRGs, Rupke, Veilleux & Sanders 2002, 2005; Martin 2005) as well as in luminous infrared galaxies (Smail et al. 2003; Swinbank et al. 2005), Lyman break galaxies (LBGs) of (Lowenthal et al. 1997; Pettini et al. 2000, 2001, 2002; Adelberger et al. 2003; Springel & Hernquist 2003) and even in some gravitationally lensed Ly emitters at (Franx et al. 1997, Frye, Broadhurst & Benítez 2002). The typical observational characteristics for galactic winds appear in most ULIRGs (e.g. Heckman et al. 2000) and most of the LBGs at redshift (e.g. Lowenthal et al. 1997). It is obvious that they are common and may play an important role in galaxy evolution. Outflows on galactic scales can remove a substantial fraction of gas from a galaxy that would otherwise collapse into cold, dense molecular clouds and form stars (Springel & Hernquist 2003; Rasera & Teyssier 2006; Stinson et al. 2006). It has been pointed out that galactic winds can effectively carry metal-enriched gas from the galactic disc through the halo and straight into the IGM, hence being the main process responsible for the observed metallicity in the IGM (e.g. Tegmark, Silk & Evrard 1993; Nath & Trentham 1997; Heckman et al. 1998; Dubois & Teyssier 2008, but see also Gnedin (1998); Silk et al. (2012)).
The exact mechanisms behind the driving forces of a typical wind are poorly understood. Observations of mainly nearby objects show that the mass outflow rate is similar to the star formation rate (SFR), and that the efficiency of conversion of supernova (SN) energy into kinetic energy of the outflow is above 10 per cent (Veilleux et al. 2005). The outflows are multiphase in nature and are limb brightened, suggesting a hollow cone structure (e.g. Veilleux et al. 2005; Sharp & Bland-Hawthorn 2010).
SNe are generally believed to play a key role in providing the vast amounts of energy that can be observed in wind-like, but also convective galactic outflows. The SN energy is thought to ’thermalise’, and generate a pressure-driven outflow similar to the stellar wind bubble (e.g. Weaver et al. 1977). A bubble produced by a single SN will expand due to its high internal overpressure, until the pressure has dropped to a value comparable to the bubble environment. Even before, the shock front produced by bubble expansion will collapse due to cooling, furthering the Rayleigh-Taylor instability, followed by a filamentary re-infall of cold, dense gas into the hot, over-pressured bubble region. By similar mixing processes, smaller bubbles are practically dissolved before they can reach the edge of the galactic disc. In order to drive a wind by pressure, it is therefore inevitable to create a sufficiently large region of thin and hot over-pressured gas, which harbours enough internal energy to provide a steady phase of expansion until the bubble can break out of the cold and dense gas disc (e.g. Mac Low et al. 1989; Stil et al. 2009; Wünsch et al. 2011). Simulations of the ISM (e.g. de Avillez & Breitschwerdt 2004, 2005) feature these overpressure regions within the hot phase of a multiphase ISM.
Recent simulations have attempted to connect the ISM structure and dynamics to mass loss from galactic discs: Hopkins, Quataert & Murray (2012) present galaxy scale simulations with parsec-scale resolution, including the effects of radiation pressure (local as well as long-range), stellar winds and supernovae. They show that different feedback processes are responsible for local and large-scale effects, leading to an ISM structured on large scales into filaments and superbubbles. Spiral structuring is also important for the ISM in some of their models. They show that the mass-loss rates from the disc depend strongly on these details within the discs, which are therefore clearly a key issue in understanding the galactic wind phenomenon. Similar to Dalla Vecchia & Schaye (2012), who study mass loss from galaxy discs with SNe-feedback only, they do not model the interaction with the gaseous halo. Both, Dalla Vecchia & Schaye (2012) and Hopkins et al. (2012) find mass outflow rates in excess of the star formation rate, in contrast to Dubois & Teyssier (2008), who do include the interaction with the gaseous halo (albeit with a simpler feedback model) and find mass loss rates of only a small fraction of the star formation rate. This supports the conclusion of Dubois & Teyssier (2008) that the interaction with the gaseous halo is another key factor in galactic wind studies.
The energy input from a star-forming region extends over a few yrs (e.g. Voss et al. 2009), which coincides with the estimated age of observed superbubbles (e.g. Bagetakos et al. 2011). Superbubbles quickly come into pressure equilibrium with their surroundings, once the energy input ceases (e.g. Krause et al. 2013). This suggests that buoyancy has a potential to become important. In elliptical galaxies, especially the central ones in clusters of galaxies, it is clear that light bubbles rise buoyantly through the hot atmosphere, dragging along cooler X-ray gas (e.g. Churazov et al. 2001; Forman et al. 2007; Roediger et al. 2007; Tremblay et al. 2012), which is enriched in metals (e.g. Heath et al. 2007; Kirkpatrick et al. 2011). For disc-dominated galaxies, which are the main focus of the present paper, the hot gas has only been observed out to about one scale height, roughly 1-2 kpc, above the disc (compare the review by Putman et al. 2012). Hence, no rising bubbles could be observed so far. From a hydrodynamics point of view it is unclear if such bubbles should be present in disc-galaxy haloes: the bubbles are likely jet-blown in the ellipticals, and would be related to star-cluster outflows in the star-formation dominated disc galaxies. Also, stability is an issue (e.g. Kaiser et al. 2005), and it is currently unclear if the expected differences in Cosmic ray content or magnetic-field dynamics (Gourgouliatos & Lyutikov 2012) would allow for stably-rising bubbles in the haloes of disc galaxies. More likely, some mass entrainment and mixing takes place (Creasey et al. 2013), reducing the superbubble entropy, such that only lower halo altitudes may be reached. ISM structure, overpressure, turbulence and possibly buoyancy might therefore be considered to be important factors for the onset and overall evolution of galactic winds.
Over recent years, a wealth of information has been collected specifically for LBGs: Their masses are of order - (e.g. Pettini et al. 2001; Steidel et al. 2010). The star formation rates are typically of order 10 yr and the star formation densities are of order 1 yr kpc (e.g. Erb et al. 2006); the latter is similar to local starburst galaxies. LBGs are well-known for their associated absorption lines (e.g. Kulas et al. 2012; Law et al. 2012), commonly interpreted as expanding (sometimes also infalling) shells of gas (e.g. Richling 2003; Krause 2005; Tapken et al. 2007; Verhamme et al. 2008; Schaerer et al. 2011). Such shells may be interpreted as radiative shock waves related to a galactic outflow. For this to occur, the cooling time of the shocked gas must be below the one of the unshocked gas, which may be arranged for by a stratified gaseous halo. Because outflow velocities are found to be typically about 150-200 km s (Verhamme et al. 2008), the sound speed in the gaseous haloes of LBGs may not exceed this value, limiting the halo temperature to K (compare Krause 2005). It should also not be much lower, because otherwise one would also expect many slower shells. In this framework, the independence of the shell velocities from the star-formation properties (Steidel et al. 2010) may be understood in terms of a wind mechanism which is just able to power a weak shock, such that the outflow velocities of the shells are always constrained to be near the sound speed in the halo. The neutral hydrogen columns of the shells are between cm and cm (Verhamme et al. 2008). Because the shells consist of swept-up and cooled gas from the galaxy’s halo, this number also indicates the column of hot gas before the passage of the shell.
Galactic winds have recently been simulated by Dubois & Teyssier (2008) and Powell et al. (2011). Dubois & Teyssier (2008) model their galaxies as cooling and collapsing Navarro-Frenk-White (NFW) spheres, and focus on the onset of a galactic wind working against the ram pressure of the in-falling halo material. They find that galactic winds arise only in low mass systems with comparatively small ram pressure which arises due to cooling and subsequent halo contraction, whereas larger ones will typically exhibit galactic fountains instead. Powell et al. (2011) study high-redshift galaxies () which are still in a phase of strong accretion by filamentary inflow of cold matter. They investigate if galactic winds may significantly alter the mass accretion rate of young galaxies in order to inhibit their further growth. They conclude that, though in these violently star-forming systems strong winds will develop, the accretion rate will not be affected, and hence there will be enough gas supply for long-lasting, intense star formation. Most recently, Verhamme et al. (2012) have performed radiative transfer on hydrodynamics simulations similar to the ones in Dubois & Teyssier (2008). They find that line profiles with enhanced red wings, similar to outflowing shells, may also be produced from features of the ISM within or close to the galactic disk, as large-scale shells are not produced in the simulations. Interestingly, their result is strongly dependent on inclination, in stark contrast to observational results by Law et al. (2012). The latter might indicate that the underlying hydrodynamic model is not representative of the majority of LBGs.
In these simulations, the simulated galaxies emerge self-consistently from some physical input. While this is extremely useful for putting LBGs into the Cosmological context, there are apparently difficulties to connect to the observations. Here, we use a fully tunable galaxy model as initial condition. We motivate the choice of our parameters from observations and Cosmological simulations and investigate the effect of the feedback strength and implementation. In particular, we make the important intermediate step to use superbubbles from star-cluster outflows instead of single supernova bubbles, and investigate the outflow properties as a function of the size of the superbubbles, thus parametrising the ISM physics which leads to the development of different superbubbles in different types of galaxies.
In Section 2 we present theoretical considerations about the most likely wind drivers, involving an analytical model to sketch the processes during buoyancy-driven bubble expansion, the multiphase ISM and the energy requirements for kinetic wind driving. Section 3 contains the simulation setup as well as the most important physics. We present our galactic wind simulations in Section 4 and compare mass and energy outflow rates for different assumptions about the stellar feedback, and the thermal pressure of the gaseous halo. We find that galactic outflows are stronger for more concentrated supernovae, less halo pressure, and if we include a thermal energy component with the supernova events. We discuss these findings in Section 5.
Here, we consider thermal energy injection, secular accumulation of kinetic energy and buoyancy in detail, and in particular their relevance for driving galactic outflows. Since the energy balance of a galactic outflow depends critically on the surrounding halo, it is necessary to set up a self-consistent model of the latter in the first place.
LBGs with winds typically occur at redshifts between 3 and 4 (compare Section 1, above). Let us therefore consider a NFW halo in hydrostatic equilibrium at redshift . The critical background density of baryons in the intergalactic medium (IGM), , can be obtained via
(Ohta, Kayo & Taruya 2003), where
with . For a flat Lambda Universe it follows from the Friedmann equations that
By choice, the model system shall have a virial radius kpc. With at , this immediately yields a scale radius kpc by invoking a value for the concentration parameter , which is verified by Zhao et al. (2009) for our underlying redshift. The baryonic mass confined within may be pinned down via the critical baryon density, , to a value of . As mentioned above, we assume that the gas in the halo is initially in hydrostatic equilibrium, and isothermal, suggesting a radially exponential distribution of baryonic matter:
where is the proton mass, and
being the NFW potential dominating at larger radii. The other two potential components are due to the disc and the central bulge, respectively, and will be explained in sect. 2.1.2. The density distribution according to equation (4) is visualised in Figure 2. Note that the density in the inner parts of the halo remains within reasonable bounds due to our simulation domain being cut off at . Since the halo shall be isothermal, we can vary thus that the density at the inner edge is not higher than typical disc density values, which are of order . The resulting temperature is K (close to the one inferred from observations compare Section 1), and by integrating the now well-defined baryonic density profile, we obtain a baryonic mass of being situated in the hot halo. If this halo gas mass would be distributed across a wind shell of 10 kpc radius, it would produce a neutral hydrogen column of cm, also what is required by the observations (compare Section 1). We have also checked that these parameters are in reasonable agreement with Cosmological simulations, which successfully reproduce observations of high-redshift galaxies (Sommer-Larsen 2006; Razoumov & Sommer-Larsen 2007; Greve & Sommer-Larsen 2008; Laursen et al. 2009; Sommer-Larsen & Toft 2010). From these simulations, we also find that the overall baryon fraction for comparable galaxies is always around or above the Cosmologically expected value, even for the strongest feedback. This is in agreement with a picture in which the so-called missing baryons are actually located in the gaseous haloes of galaxies (Sommer-Larsen 2006; Ntormousi & Sommer-Larsen 2010). For the Milky Way this picture has now been confirmed observationally (Gupta et al. 2012). In the Cosmological simulations, we find about 10-20 per cent of the baryons in the hot halo. We have chosen here a smaller value of 5 per cent, still consistent with observations, which should favour the development of winds. This would also account for some halo gas being lost already due to preceding galactic wind activity. We also note that the stellar to halo mass ratio in our simulations, 7.3 per cent, is well within the range given by Moster et al. (2013) for statistical halo abundance matching for our halo mass and redshift, 0.08 - 25 per cent, but towards the high side of their central value, 1.4 per cent. We discuss the implications of these choices in Section 5.
The larger part of , still amounting to , must be considered to have settled into the disc. With the halo density given for all radii, the halo pressure can be obtained from the ideal gas equation
where , due to ionisation. The initial equilibrium state for the halo will only hold as long as the temperature is kept constant. Yet since in some of our runs radiative cooling is permitted for the model halo, the subsequent temperature decrease leads to some contraction of the halo with time. This in some sense accommodates for the fact that galaxies at the given redshift are still accreting halo material in significant amounts. However, the interaction between (filamentary) infall of material into a galactic disc and the onsetting wind is beyond the scope of this work and is studied thoroughly by Powell et al. (2011).
The task of constructing an isothermal halo in hydrostatic equilibrium is encumbered by the condition that its density should converge against a certain background value. A halo potential of the form
given by a constant rotational velocity for large may seem physically justified. Such a model is described closer in Flynn, Sommer-Larsen & Christensen (1996), however, this model entails the fact that the halo pressure will not converge. This means first of all that shock fronts could theoretically proceed to infinity as due to the resistant pressure decreasing strongly with they will accelerate forever. Furthermore, the density would have to drop adequately in order to maintain a constant temperature all over the halo, and would soon reach unreasonable values below the cosmic background (compare Figure 1). We hence adopted the NFW potential to overcome the described problems.
Several approaches to establish a stable disc-halo system have been tested. A detailed description for a possible setup can be found in Cooper et al. (2008). In general, the following issues have to be kept in mind: Firstly, we want the gaseous disc to be rotationally supported (i.e. in hydrodynamic equilibrium), whereas the halo shall be pressure-supported (i.e. in hydrostatic equilibrium), which inevitably causes friction and shear effects in the transition zone. In addition, the halo cannot be truly set up in a pressure equilibrium with the disc, as the halo isobars are geometrically not parallel to those of the disc, which inevitably causes some motion in the halo. Therefore, we allow the system to relax for one Myr. The resulting setup is then sufficiently close to an equilibrium configuration to allow for the development of relatively stationary outflow solutions (compare below). As mentioned above in equation (4), the total potential is built up of three components, where the disc component is a combined form of a Miyamoto-Nagai potential (Miyamoto & Nagai 1975):
The bulge component is basically a central potential,
These two components are further described in Flynn et al. (1996), which we will use as the basic prescription for our disc setup. We have scaled down the mass-related parameters therein (, , , and ) by a factor of 0.37 to match the residual disc mass (gas and stars) of . The length-related sizes (, , , , and ) in the description by Flynn et al. (1996) have been scaled down by a factor of 0.5 for our purpose, leaving our disc at a scale radius of . An overview of all the values related to the bulge and disc potentials is given in Table LABEL:tab-discflynn. For comparison, a typical LBG is observed to have comparatively small size, and a mass probably an order of magnitude smaller (a few ; see Pettini et al. (2001)) than the more massive SINS galaxies (Genzel et al. 2008). They form stars dominantly in a steady mode with a range of star formation rates, tens of solar masses per year not being uncommon (Pettini et al. 2001; Springel & Hernquist 2003). With a gas fraction of 50 per cent in the disc, the star formation density ranges from 0.06 yr kpc to 1.4 yr kpc at the inner edge of the disc. The average value is 0.3 yr kpc, towards the lower end of the values reported by Erb et al. (2006).
Having 50 per cent of the disc mass locked in stars gives us some freedom of choice for the gas density distribution, because the disc potential is made up by the combined mass of gas and stars. We use an exponential (in radius) gas density profile with a cutoff at , that is vertically non-stratified. This latter fact is unproblematic since the disc will be given enough time to relax, so stratification will develop in the early course of the respective models (). The disc density thus reads
with being the disc scale radius, and . The disc has a vertical height of 500 pc, and thus the total gas mass is , i.e. about 50 per cent of the mass implied by the disc potential. As an alternative to the exponential distribution, one could use a constant gas density profile, which has been observed e.g. by Bendo et al. (2010) for NGC 2403; this is to be dealt with in a future paper. The disc gas pressure follows from the ideal gas equation (7), just as for the halo gas pressure.
The gravitational force is accounted for by the implementation of as an external potential. In Figure 2, the resulting density for our disc-halo system is shown as a contour plot; the initial and boundary conditions will be further explained in sect. 3.3. This setup condition applies to the complete set of simulations presented here and listed in Table LABEL:tab-sims.
2.2 Overpressured superbubbles and the multi-phase interstellar medium
Single SN bubbles dissolve into the surrounding ISM on a scale of order of
10 pc, which may also be derived from the LBG simulations presented below under
the assumption of pressure equilibrium. This is a small value compared to the scale
height of disc galaxies (compare Section 1). Superbubbles from clusters
of stars formed by the winds and SN explosions of many massive stars are therefore
more promising drivers of galactic winds. Superbubbles with energy requirements of
tens to hundreds of SN and diameters of hundreds of parsecs are frequently found
in nearby galaxies (e.g. Bagetakos et al. 2011, also Section 1).
Using such observed superbubbles directly in simulations also avoids the
possibly complex interaction of smaller interstellar bubbles leading
to the emergence of the superbubbles in the first place (e.g. Krause et
We adopt this approach here. With our recipe for superbubble injection explained in Section 3.2, superbubbles with prescribed energetics are injected on a scale of about 100 pc. This leads to a dynamic, multiphase ISM, structured on large scales with hot bubbles and cold filaments.
2.3 Global kinetic energy
A galactic wind could also be launched by the sheer amount of kinetic energy which accumulates within a gas-rich galactic disc over time. Let us consider a galactic gas disc with a mass of the order , as commonly found for Lyman-break systems. Let us further assume the SFR to be around , which means in turn that we are going to encounter about one SN every 10 years. Normalised to the entire mass of the system this would mean a SN rate 10 times as high as in the Milky Way, which too has a gas mass of order (the larger part of its mass is locked in stars) at a SFR of . Our model system may therefore be regarded to be in a starburst phase. SNe are known to give rise to considerable turbulent motions within a disc (Dib, Bell & Burkert 2006), each yielding a contribution of erg at a presumed efficiency to the overall kinetic energy stored within the gas phase of its host galaxy. Unlike internal energy, kinetic energy has the advantage that substantial fractions will not be radiated away immediately, but rather dissipate on the dynamical timescale (Mac Low et al. 1998; Burkert 2006). Allowing the turbulent energy to pile up for would result in an energy reservoir of order erg for the disc as a whole. Since the gravitational binding energy is known to be of the same order for a system of kpc radial extent, material ejections from the disc into its surrounding galactic halo indeed becomes plausible at a certain point in time. The approach of launching a turbulence-driven outflow has been investigated by Scannapieco & Brüggen (2010). In their models SN feedback was simulated by injecting unresolved kinetic energy, which is described by an isotropic pressure term in the Euler equation. Here, we also investigate kinetic energy driving, studying models where we inject only kinetic energy, instead of a combination of thermal and kinetic energy (compare below). We will however resolve the kinetic energy.
2.4 Buoyancy of superbubbles
Even large superbubbles may come into pressure equilibrium with their surroundings while still being close to or even within their galactic disc. In this case, buoyancy needs to be considered. An interesting physical quantity in this context is the entropy index . Here, we calculate at relevant locations within the underlying NFW halo at redshift . The entropy index is defined as
where is the number density of particles, the units of being given in keV cm. Generally, a bubble with an entropy index higher than its environment will experience a buoyant force, meaning that with S being known everywhere, we can easily determine the height a buoyant bubble can reach.
With the pressure expression from eq. (7), the entropy index transforms into
Let us consider a bubble produced by a single SN in an early state of evolution. The entropy index is highest within the central hot gas phase of the bubble, and this is the region most relevant regarding buoyancy. Note that S is defined such that during the process of adiabatic expansion it is not going to change over time. For the hot bubble interior, S may decrease due to mixing and cooling. Cooling times are long compared to the simulation time, and mixing shall be neglected here in the first instance. This in turn means that the phase of evolution in which we investigate a bubble doesn’t matter all too much. A typical SN will release about of energy. From the equation of motion for a blast wave in the thin shell approximation, it follows that 60 per cent of this energy will be in the form of thermal energy. Implying an ejecta mass of 8 and a bubble in an advanced state, e.g. with a radius of to start with, the density will be of order . It follows then, assuming a temperature of , that the entropy index from eq. (13) reaches several . Given a typical entropy index for the gas disc of order , the former value is certainly enough to raise the bubble away from the disc midplane into the disc-halo transition region. In our example, the values for and in the halo amount to and , respectively. Hence inside a bubble formed by several SNe will be typically high enough to exhibit buoyancy effects within the halo at least at low radii. This conclusion might however be affected by the (unknown) mixing of the different ISM phases. In our simulations, we include the buoyancy effect of the superbubbles. We inject the bubbles with even higher entropy index, because numerical mixing - we have to inject the superbubble on a scale of a few grid cells - strongly reduces the entropy index. The energetic effect of buoyancy is however likely minor: For a sized bubble and typical parameters of our simulations below, the acquired velocity would only be about 100 km s. Thus, the expected effect is that superbubbles first expand to pressure equilibrium and then hover near the disc-halo interface. If fed sufficiently by other superbubbles, they may develop into a galactic wind. Otherwise, the bubble overshoots, collapses again and dissolves.
3 Numerical Methods
We perform 3D simulations with the magnetohydrodynamics code NIRVANA (Ziegler & Yorke 1997) on a spherical grid.
Thus angular momentum is well conserved. In most runs, we simulate only a small fraction of the azimuthal angle, which allows to explore a larger part of the parameter space.
We have parallelised the code making use of the Message Passing Interface (MPI) library. Our simulations run for typically 48 hours on 6 SGI Altix processors. The gas evolution is calculated by solving the continuity, energy and Euler equations. A constant background gravitational potential accounts for the stellar and gaseous disc components, a bulge and the dark matter halo (compare above). The radiative cooling function used here is the equilibrium cooling curve described by Sutherland & Dopita (1993). It accounts for the overall metallicity which is assumed to be equal to the solar metallicity, and operates only within a temperature range between a lower limit of K and an upper limit of about K, with the exact value depending on the respective halo equilibrium temperature: For some of our models, the only effect of the upper cutoff is to prevent cooling in freshly injected SN shells. This is required to establish a more realistic SN remnant, before the shell cools and the remnant enters the snow-plough phase (see Sections 2.3 and 3.1 for more details). We also use the upper cutoff to entirely inhibit cooling of the halo in some simulations. We do this to account for the unknown halo metallicity, which has a strong impact on radiative cooling. In this way, we cover the limiting cases of strong and negligible cooling of the halo.
We shall briefly describe the most important methods used in our studies and, if non-trivial, justify them physically; this includes cooling restrictions, SN triggering and their blastwave implementation.
3.1 Stellar feedback
Star formation is triggered randomly for each single cell as soon as certain criteria are met, and SNe occur immediately in an amount related to the mass of stars produced. Star formation criteria include a local surface density exceeding the critical value pc required for star formation to set in (Kennicutt 1998). Before calculating the local surface density, a volume density criterion applies for each cell to ensure that it is part of a region dense enough to produce stars, which is, in particular, the disc. Cells having a density less than are considered to be either halo cells or too rarefied for star formation to set in. In a few special cases, large high-density gas regions can be found far away from the disc. We are to assume then that our model galaxy is essentially breaking up as a consequence of too strong feedback. In consequence, once the disc has lost integrity, the Kennicutt-Schmidt law might no longer apply. To avoid perturbations from this effect, the column from which the surface density is calculated, comprising only the aforementioned disc cells, shall be no higher than one fourth of the total -range of the simulation domain. This value chosen here, however, is not a critical parameter. Finally, in order to allow the system some relaxation after setup, SNe shall not occur before .
In this model, we regard only SNe type II, since star forming galaxies are observationally dominated by this type. Given the Salpeter IMF for the stellar mass distribution, we can easily calculate that of of gas locked up in stars, one type II SN progenitor exists, with the latter typically being as massive as on average, considering stars within a range from 8 to 120 .
We assume that stars in all our model galaxies form in accordance to a local Kennicutt-Schmidt law (Kennicutt 1998), given by
where denotes the respective surface densities for star formation and gas mass. The gas surface density is calculated for every time step and every grid point within the --plane by integration of all disc cell masses along and dividing by the surface area of the respective column. Integrating along instead of the normal with respect to the disc midplane is a sufficient approximation since the disc extends only across a small angle . Moreover, constraining star formation by limiting the maximum column height to one fourth of the -range of the simulation domain will ensure that the angle of integration is sufficiently small.
The SNe in our simulations are assumed to cluster in groups of 20 to 200 at a given time and place (henceforth referred to as a ’SN event’). Thus, we can calculate a probability value for every disc cell and each time step, giving the likelihood for a SN event comprising SNe,
with being the number of disc cells in the respective range of integration along . , referred to as the ’event size’ herein, is a preset parameter which will be kept constant during each single simulation. A random number is then drawn for each disc cell at every time step. The occurrence of a SN event is then triggered according to the local probability. Note that any altering of the resolution has to come with an appropriate change in the event size range; too high event sizes will require a manifold of the gas mass available in the cell, too small event sizes may produce unresolved bubbles.
3.2 Blast wave implementation
In this subsection we provide details of our blast wave implementation and follow the evolution of a single superbubble in a test simulation. If a SN event is determined to occur for a specific cell, the following modifications in mass and energy will immediately take place: An amount of is regarded to be no longer available in gaseous form since it is bound in stars, and henceforth removed from the cell. This amount can, in case of large , exceed the cell mass, however, within our range of event sizes the total mass deficit due to this error is below a five per cent threshold for the entire disc, and hence considered negligible as it will insignificantly alter the disc dynamics. 25 per cent of the newly formed stellar mass is returned to the gas phase due to stellar winds and SN ejecta. So, essentially, our code removes of gas from the SN-triggering cell. The remaining mass is distributed equally among the six neighbouring cells except for a small remainder of g/cm within the central cell, so that the density increase is the same in all six adjacent cells. We assume an energy injection of per SN. According to the thin shell approximation for blastwaves with instantaneous energy injection, 60 per cent of this energy is released as internal energy, fed into the SN cell and thus building up an overpressure with respect to the surroundings. The remaining 40 per cent of the energy total is kinetic energy, added as an extra velocity component to the neighbour cells. This velocity of the SN ejecta, , is typically greater than upon release, and therefore supersonic with respect to the sound speed inside the dense disc material, in agreement with superbubble observations (e.g. Bagetakos et al. 2011). The temperature in the SN cells varies due to the spherical geometry of the grid. It is however always of order K or higher, exceeding the value determined by Dalla Vecchia & Schaye (2012) to achieve a converged outflow behaviour, . While this high temperature prevents cooling in the bubble interior, cooling is still very efficient in the shocked shell of the injected superbubble: In an interstellar bubble, the energy lost by adiabatic expansion from the bubble interior is used to accelerate the surrounding shell, which drives a shock into the ambient medium, where the energy finally thermalises and radiates. The total energy of a superbubble is lost within about yrs, once the energy input has stopped (Krause et al. 2013). Thus, in order to establish superbubbles of hundreds of parsecs diameter on the grid, as observed, it is necessary to suppress cooling in the shocked superbubble shell for a short time interval after injection of the superbubble. This is not unphysical, since our simulations are not meant to explain the origin of the ISM structure, for which radiation pressure and stellar winds are crucial (Hopkins et al. 2012), but to explore the effects of superbubbles with given sizes. We implement this by a threshold value above which no cooling is taking place. For most simulations, we have chosen this threshold value to be K, slightly above the halo temperature. In this case, no cells apart from the shells of freshly injected superbubbles are affected. Only in simulations with non-cooling gas haloes we artificially suppress cooling of the halo gas by setting the threshold to a value slightly below the halo temperature, in order to also suppress cooling in the halo gas. This initialises slightly bigger superbubbles in the latter runs. But the effect is shown below not to be significant, as the mass and energy loss rates are higher for the case when the halo is allowed to cool.
The choice of K for the threshold value implies a critical shell-expansion velocity of 271 km/s, above which the expanding shells of freshly injected bubbles are assumed here not to cool. Hence, cooling is allowed way before a superbubble reaches observed velocities ( km s, Bagetakos et al. 2011). A much higher threshold value (say factor of ten) would push the critical shell velocity up, well into the regime of single-supernova shells, which are not addressed in our simulations. If applied to our superbubble setup anyway, the injected energy would be radiated away quickly, and the bubbles would dissolve without producing any large-scale effects. Because for blastwaves the postshock temperature depends on the radius to the third power, smaller ( factor of two) changes in the threshold temperature do not change the physics significantly. We have determined the threshold value experimentally to ensure that the superbubbles are just established properly on the grid. The exact choice of the threshold value therefore has an effect on the net energy input. However, we do not model the emergence and early evolution of the superbubbles (compare e.g. Krause et al. 2013, for a discussion), but use idealised superbubbles instead. We study the effect of variations of the superbubble-injection mechanism in Section 4.2. However, instead of varying the threshold temperature, we directly vary the amount of injected energy per bubble.
To check the behaviour of our blast wave implementation, we have modelled a box of cells on a spherical grid section of , and in each and direction. There is neither an external potential, nor does any other force (e.g. centrifugal) apply. The overall density is set to , and the temperature to which is a common value for disk material in the models presented below. An energy equivalent of 100 SNe, or erg, is released at right in the centre of the box as described above, forming an over-pressured, expanding hot gas bubble within a few 10,000 years (Figure 4). There is no cooling taking place in this test run. We follow the bubble expansion over , tracing the distance between the shock front and the centre of explosion (Figure 3) as well as the energy decrease with time. The shock front in the underlying model is found to expand with a law, as expected. Note however, that this is the expansion behaviour as expected from a bubble produced by one single SN. Superbubbles powered by many SNe spread out in time should rather expand with (Oey 2009). This is because all of our bubble-producing SNe are triggered in one cell within one time step, as resolution prevents us from spreading SNe reasonably in space and time, in order to produce more realistic superbubbles. We estimate that this effect increases our bubble sizes artificially by about 25 per cent despite the smaller expansion rate, as we start out with a much higher energy. On the other hand, we find that about 10 per cent of the initially released energy is lost by numerical effects within the first 100,000 years, however, any further loss thereafter is comparatively small. Because the advection step of our code conserves only the thermal and not the kinetic energy exactly, preferably the kinetic energy component will be lost. Figure 4 shows two snapshots of the SN bubble evolution, respectively 0.1 and 2 Myr after the event was triggered; the inner, rarefied region carries the internal energy. The kinetic energy resides within the compressed high-density region surrounding the bubble. Its slightly asymmetric form and imbalances in the kinetic/thermal energy distribution are a result of the coarse implementation.
3.3 Setup, boundary and initial conditions
|Resolution111For the angular coordinates, the resolution at the inner and outer radial boundary is given.||SN energy222The energy released by a SN event is subdivided into a kinetic and a thermal component.||Event size||Res.333The term ’Residual density’ refers to the density left over in a cell after being subject to a SN event.||Halo|
We run our simulations on a 3D spherical grid, with the radial dimension extending from 0.4 to 10.2 kpc, the polar angle covering a section between and , and the azimuthal angle covering only a narrow ’wedge’ of the disc within and in range. Note that the space close to and as well as the one at must be omitted, as due to the spherical geometry grid cells within this space would become increasingly narrow. This in turn would lower their crossing timescales significantly, requiring high computing times for the innermost zones. The simulation of just a small azimuthal sector of the disc instead of the whole -range implies large-scale rotational symmetry. We have also performed control runs relaxing the assumptions about the azimuthal extent (STaz) and the polar-axis cutout (STpol). The polar axis cutout affects the energy fluxes by about 10 per cent. Otherwise, the effects are minor and are discussed in the appendix.
For the main runs, the simulation domain is divided into grid cells in , and directions, respectively. Thus, a region near the disc midplane at is spatially resolved to
(compare Table LABEL:tab-sims for details). Our choice made here concerning the resolution will be explained in more detail at the end of this section. We choose reflective boundary conditions for the lower boundary and the upper and lower boundaries each, whereas on the outer boundary in -direction inflow and outflow of material shall be permitted. The boundary conditions for the boundaries in azimuthal direction () are chosen to be periodical. Table LABEL:tab-sims shows the total set of simulations performed in the frame of this work with their respective parameters.
We have investigated the resolution dependence of the SFR (see Section 5.1 below for a discussion of the dependence of the outflow rates on resolution), varying the reference resolution at radius from 16 pc to 65 pc (R16 - R65, compare Table LABEL:tab-sims) for a standard simulation. Assuming one SN in 100 of stars formed, we find the SFR in our system at all resolutions to be about , yielding an SFR per unit mass of . As a comparison, this is several ten times the SFR per unit mass in the Milky Way, which would be of a few . Our SFR is therefore in the relevant range; e.g. Pettini et al. (2001) observe values of about for their sample of -LBGs at redshift , which, accordingly, would result in an SFR several per unit mass (or a few SNe per year). The overall SN rates of our model galaxy are displayed in Figure 5 for all resolutions. The graph for resolution shows the strongest deviation, indicating that too coarse resolutions will notably affect the star formation rate. All graphs agree within 26 per cent, however, if we regard only resolutions of and finer, the error reduces to nine per cent.
We begin with an investigation of how the method of SN energy injection affects the emerging wind. For this purpose, we have run a set of simulations with . One simulation uses the Sedov-Taylor blast wave model, and hence both kinetic and thermal energy are injected with every SN event (denoted ’ST100’). In addition, two models were calculated, injecting a purely thermal energy fraction of 40 per cent (denoted ’TE0.4’), and 60 per cent (’TE0.6’) of the total SN energy yield, respectively, and another one, injecting a purely kinetic energy fraction of 40 per cent (’KE0.4’). The characteristics of the pressure-driven and the kinetic energy driven cases are discussed in the first two subsections, respectively.
All the runs presented in subsections 4.1 and 4.2 include a cooling halo. Since halo pressure is reduced by cooling, winds will arise comparatively easily in this case, allowing for more prominent effects more suitable for later comparison. Subsection 4.3 investigates the question how the sizes of SN bubbles can affect the strength of galactic winds; for this we have run another set of three simulations featuring Sedov-Taylor blast wave models and different event sizes each. In contrast to the previous runs, the runs in subsection 4.3 are each performed twice, with both, a cooling and a non-cooling halo, respectively, to investigate the limiting cases of the possible effects of varying metallicities in such objects. We show that the different halo pressures have a significant effect on the wind. All of our results herein will then be compared in the final subsection.
4.1 Pressure-driven winds
In Figure 6 we show the mass density distribution of simulation ST100 at times of 10, 20, 40, 60, 120 and . We can clearly discern individual superbubbles expanding already at beyond a height of above and below the disc.
At this point they have expanded out to pressure equilibrium and hover above and below the disc due to
the long buoyant rise time. Bubbles in the outer part of the disc collapse back to the disc.
Bubbles inside a radius of 3 kpc however merge and are fed sufficiently
quickly to prevent them falling back to the disc. The region of the disc where
this happens has a star formation density greater than about 0.1 yr kpc.
These bubbles keep expanding, driven by their overpressure against the radially quickly declining halo pressure. At the superbubbles
create a low density funnel close to the axis of symmetry. Since the gas inside this structure provides less resistance to subsequently escaping superbubbles than the rest of the halo region, material from the succeeding bubbles will continue to flow at ease through the funnel.
The individual expanding superbubbles are however still identifiable by individual shock fronts, which may be more easily seen in the accompanying movie of run STaz. The funnel
is surrounded by a conical structure of notably denser material which was originally entrained from the dense disc by outgoing shock fronts and hence continues to move outwards. Over time, enormous amounts of SN energy are fed into the disc, which in turn becomes extremely turbulent: large portions of gas are torn out of the disc midplane, partially due to entrainment by the wind, but eventually fall back. The shape of the disc gets highly irregular and clumpy but the disc remains overall intact.
Since we are dealing with a rather massive system, it might seem likely, regarding the studies by Dubois & Teyssier (2008), that outflows appear preferably in the form of galactic fountains. We make here the usual distinction (compare e.g. Dubois & Teyssier 2008) between the two common types of outflow solutions: A wind is defined to be supersonic with respect to its internal sound speed. A fountain, on the other hand, is subsonic. Galactic fountains are therefore much more susceptible to the Kelvin-Helmholtz instability and usually turbulent. Both types of solutions may in principal be bound to the galaxy or reach escape velocity. The smaller bulk velocity of the fountain solution usually prevents it from escaping the galaxy and the flow becomes convective, lead by a roughly spherical weak shock or sound wave around the whole system.
In contrast, the bulk velocities in the wind gas may easily reach escape velocity. Due to the geometrical constraint from the galactic gas disc, the outflow becomes conical. Figures 6 and 7 demonstrate that the outflow which has emerged in run ST100 has developed all the usual characteristics for a wind solution. The escape velocity at 10 kpc distance from the disc amounts to , which is well below the typical wind velocities close to . The difference to Dubois & Teyssier (2008) is mainly the size of the disc. Dubois & Teyssier (2008) have chosen a much larger disc and therefore might not reach the required SN density to drive the outflow.
4.1.1 Mass outflow
For a quantitative analysis of our models, we calculate the net mass flux across a spherical shell of inner radius and outer radius . We start with
which is the net mass flux at any point in time for a spherical layer of grid cells at a given radius . The factor is a correction term which accounts for the fact that our simulation box covers only of the total range. Due to the box limits in range, a part of the wind at the poles is neglected. Due to the small surface area, this error is not significant (of order 1 per cent, compare appendix). The average mass flux for all layers at radii is determined every 1 Myr, and then averaged over 10 Myr, yielding the total net mass flux . Figure 8 shows mass flux rates from 0-200 Myr for run ST100 across shells of respective thickness of for various shell positions. In the innermost shells, winds show up earlier and stronger, however, a large fraction of the outflowing mass in these inner shells is likely to represent entrained disc material. This material might, in some cases, fall back soon after its ejection from the disc, and actually not contribute to the mass carried away by the wind.
4.1.2 Energy outflow
To obtain the net energy flux, we assume the same shells as before. The energy flux comprises kinetic and thermal components:
The mean value for the net energy flux is averaged in the same way as the net mass flux, namely . Here, the polar contribution is a bit higher and we show in the appendix that we underestimate the energy fluxes by about 10 per cent due to the polar cutout. Again, the energy flux rates displayed in Figure 9 represent different shells of thickness each, for different shell positions.
Comparing the respective shells of measurement in Figs. 8 and 9, we can clearly see a convergence of the graphs with increasing shell radius. Measurements closer than exhibit more pronounced extrema, and, in case of strong turbulent feedback or irregularities in the disc, may be prone to notable perturbations arising from the disc. If too close to the box boundary at , interactions with the boundary itself might distort the actual result in a few cases. Therefore, we choose the range in between as the most reliable one.
All plots exhibit one more or less strong peak, which is the first shock front clearing the path for the wind yet to come. Any further peaks are a result of local and temporal concentrations of SN events; yet these anomalies will be mitigated as the energy outflow will stabilise over time. The basic level of energy carried by the wind is several . So, with an average input of some in our models, we can define a wind efficiency as the ratio of wind energy to injected energy. The latter is stable on a level around , as is shown in the lower panel in Figure 9.
4.2 Kinetic energy-driven outflows
In order to compare directly the respective importance of the thermal and kinetic forms of energy injection, we have performed three simulations, where we inject only thermal energy or only kinetic energy (Figure 10). Note that these simulations permit cooling in the halo, which subsequently reduces the environment pressure the wind has to overcome. The cooling halo is particularly necessary for the sake of the comparison in this section; without it an outflow may not be strong enough to leave the disc at all in some of the presented cases. In run TE0.6, we inject the thermal energy component, only, using the standard fraction of per injected SN. This run has a slightly slower wind start, but later on is statistically indistinguishable from run ST100 regarding mass and energy outflow rates (Figures 11 and 12). Using only the 40 per cent kinetic energy (KE0.4), the outflow is much weaker: It has now a much harder time to get out of the disc. The part in the hemisphere with negative values is even dragged back by the ram pressure of the infalling halo (120 Myr). The outflow stalls completely between 110 and 120 Myr (compare Figs. 11 and 12). These results seem to indicate that the thermal energy part is the more important one for wind driving. We have also performed a run (TE0.4) with the thermal energy injection being reduced to the level of KE0.4. Here, the wind is also noticeably weaker, and the downwards going bubble also comes back. The statistics indicate a stronger outflow for TE0.4. However, the system is evidently just around the threshold, where it can drive a wind at all. Therefore, small changes might affect the result strongly. Remembering that our numerical scheme conserves the thermal energy better than the kinetic one (compare section 3.2), we conclude that the differences between KE0.4 and TE0.4 are not significant.
4.3 Bubble size
The last set of simulations presented in this study features a variation of the event size introduced in Section 3.1, above. The event size specifies the number of SNe comprised in one single bubble. On average, for of newly formed stars we expect one SN and a gas mass return through stellar winds and SN ejecta of . This in turn requires a minimum available mass of per cell for . However, there is a chance for a mass deficit to occur, typically in the outmost parts of the disc where the defined minimum density of is just reached, or in cells close to the inner radial boundary which exhibit small absolute angular diameters. On the other hand, the average cell mass will be , which is well above the requirement for a 200-SN event. The mass deficit is not a severe issue, since in reality, the mass would come from neighbouring cells, and because the global error on the mass budget is small, no significant effect on the dynamics is expected. Locally, one might expect that we might artificially somewhat damp the kinematics in the large bubble simulations because of the slightly higher inertia in these runs. Yet, as we show below, we find that large bubble simulations exhibit the strongest winds.
Note that in some of the following simulations (NC20, NC50, NC100 and NC200) the threshold above which we inhibit radiative cooling is reduced below the halo equilibrium temperature of 600,000 . We include these non-cooling simulations in addition to the ones with cooling at solar metallicity, in order to investigate possible effects of metallicity: For metal poor gas haloes, the cooling time is prolonged. Such galaxies will therefore likely have a hydrostatic halo as we describe it. For increasing metallicity, the thermal pressure will drop due to cooling but at the same time ram pressure due to the inflowing gas will increase (compare Dubois & Teyssier 2008). With the approximations of solar metallicity cooling (ST) and non-cooling (NC) haloes, we try to capture the extreme cases, keeping in mind that a full parameter study in a Cosmological setup is clearly beyond the scope of this work. The values chosen for in these simulations are 20, 50, 100 and 200 SNe, respectively (compare Table LABEL:tab-sims). In the following, the total SFR, the onset of the wind and its temporal development will be of particular interest. We will further investigate the mass and energy efficiencies in the same manner as above. It may seem reasonable to assume that, since smaller bubbles are situated much closer to each other than large ones, dense material in between will be further compressed until star formation sets in, thus providing a positive feedback to the SFR. Yet, large bubbles may proof more powerful when it comes to triggering the wind, and thus we could find that a larger , though providing little less energy input, results in a slightly more efficient wind.
4.3.1 Star formation
A look at Figure 13 immediately reveals that the cumulative SFR for different bubble sizes undergoes little change within 8 per cent, until just before . This difference grows, being already around 24 per cent at . An explanation for this could be that large bubbles result in a violent blow-away of large gas portions, whereas small bubbles, due to their numerous occurrence, smear out the disc material over a comparatively large volume, reducing the chances for the gas to pile up in high amounts on any single spot. Both effects can result in a visible reduction of star formation, and hence the optimum range for star formation comes to lie in between 50 and 100 SNe per event.
In Figure 14 we plotted the mass-weighted height of the gas above the disc midplane, which calculates as
The resulting value for indicates the average height of all gas portions in above the disc plane at any given time. We find that for NC200 is significantly larger than for NC100, but only between and , while NC20 and NC50 show comparatively little difference. NC100 and NC20 however increase strongly in the last . Increasing values mean that during this time much of the gas is torn out of the disc forming filaments, which constitute large quantities of gas unavailable for star formation. But if this were to be the reason for the lower SFR in NC200, we would expect the NC200 graph to dominate clearly from about onwards. This possibility can hence be excluded.
In contrast, small bubbles of 20 SNe should have a smoothing effect on the overall density profile of the disc. The number of columns with respect to their density is visualised in Figure 15, whereas the total column number includes all columns within and is integrated over the total simulated time span of . The curve for NC20 should exhibit more moderate values than its NC200 counterpart, whereas extreme values below and above about should be less present in the former. Columns of high density contribute most of all to the global SFR, and should be most present in the NC50 and NC100 curves. We find however, that neither of the four curves matches any of the expectations. Therefore, we can also exclude smoothing effects inside the disc from large numbers of small bubbles to be of notable effect to the SFR. This means that the SFRs in our simulations are set by a more complex interplay of processes.
4.3.2 Mass and energy flux
Figs. 16 and 17 show the absolute and relative mass flux, and the absolute and relative energy flux, respectively, for the different bubble sizes. It has to be borne in mind, that after the differences in the SFR become somewhat stronger (compare section 5.1 below). There is no doubt that the mass flow curve for NC200 starts earliest, and much higher than the others. Early starting curves are a clear indicator that the wind developed fast; in the case of NC200 it takes for the wind to reach the radius of measurement at , giving it an average speed of nearly . Curves starting late suggest that the wind is setting in at a later point in time, but could also indicate a slower wind. The former case however applies to our simulations. The wind in NC100 starts early and still carries comparatively large mass. Of the two remaining ones, NC50 exhibits a stronger wind at a late start, whereas NC20 starts with little mass at an earlier time. When looking at Figure 17, it becomes more obvious, that large SN bubbles show a tendency to start blowing a wind in a powerful way. The NC200 and ST100 energy curves stay roughly constant in time, whereas NC20 exhibits a more chaotic behaviour after the onset of the wind. While the order is not strictly maintained throughout the simulation time, there is a clear general trend for larger bubbles to produce higher mass outflow rates. This is also generally confirmed from the cumulative numbers (Table LABEL:tab-flux): The two large-superbubble simulations have a mass outflow rate which exceeds the one of the two small-superbubble simulations by about an order of magnitude. Run NC200 has formed 11 per cent less stars than run NC100, and still ejects 34 per cent more mass. Only for run NC 50, we find a 21 per cent smaller outflow rate in comparison to NC20, while the star formation rate is 24 per cent higher.
The trend is even more evident in the cumulative energy outflow rate (also in Table LABEL:tab-flux): For all the NC simulations, they increase monotonically with superbubble size, even if normalised to the star formation rate.
4.3.3 Halo pressure
|Run||Cumulative mass flux||Cumulative energy flux|
In Table LABEL:tab-flux the complete set of runs ST20, ST50, ST100 and ST200 is compared to their respective NC counterparts. The displayed values are the cumulative mass and energy flux rates until , in absolute numbers, and , respectively. For each bubble size, the flux value of the respective ST run is normalised by the value for the respective NC run. It is obvious that for the smaller bubble sizes, and , the outflow is stronger in the absence of thermal halo pressure. Moreover, ST20 and ST50 feature one major outburst each, where a massive local concentration of feedback energy leads to the ejection of a large share of hot gas from the disc. However, if , a steadily blowing wind arises also for the thermally pressurised halo; we find both mass and energy outflow rates for ST200 and NC200 to range in the same order of magnitude, respectively. represents a special case, where an exceptionally large filament is torn out of the disc after , which accounts for the bulk of mass and energy (also compare Figure 18). If this phenomenon is neglected, the flux values for ST100 and NC 100 would be of comparable magnitude.
Figure 18 shows one snapshot from all eight runs at the same time of . We find for that in both cases the small bubble size only triggers a weak wind. In NC20, filaments bordering the upper and lower wind conus are absent, indicating that the halo pressure has already begun to force the wind conus back into the disc. In ST20 we find the wind to be asymmetric, being at least stable on one side of the disc. The same applies to ST50, where the wind is also dominant on one disc side only. NC50 in contrast developed a biconically stable wind, however, the conus is already in the process of being crushed. The wind in NC100 has ceased entirely; instead we can see the disc being just a few Myr before complete disruption - which also explains the enormous mass and energy outflow rates towards the end of NC100. ST100 on the other hand exhibits a clear biconical wind structure, with a wind steadily blowing in both directions. Stable winds also occur in NC200 and ST200. This supports our assumption that large superbubbles generally seem to boost the overall strength and steadiness of the wind. Furthermore it appears that for smaller bubbles the environment pressure becomes important: If the halo is thermally pressurised, winds arise but cannot overcome the halo pressure in the long term. In case of a cool, less pressurised halo, winds are on the brink of developing towards a stable, steady state; asymmetric developments with at least one of two coni being stable are not unlikely.
In summary, mass and energy outflow rates in the NC runs consistently show the same trend: if the event size is varied, the outflow rates for
large will tend to start comparatively high, and change barely over time. Small will cause the wind to set in less forcefully, and, as is the case with NC20, undergo occasional drops in strength. The efficiency of mass ejection in our NC models is typically between and unity. The efficiency of feedback energy conversion exhibits a convergence for most runs against , while values of are still common, and is already rare.
For the ST runs (low halo pressure), no clear trend can be discerned. A high halo pressure efficiently pushes smaller bubbles back into the disc, but low halo pressure enables bubbles of all event sizes to enter the halo overpressured and keep expanding. Therefore, the outflow properties do not depend systematically on the bubble size in the latter case. Instead, they tend to be dominated by single events, like the high concentration of SN bubbles leading to a violent ejection of large gas masses in ST20 and ST50.
5.1 Methodical accuracy
We have produced successful models of LBGs launching galactic outflows, in order to shed some light onto the exact mechanisms responsible for the onset of a galactic wind. These models feature a realistic galaxy setup with superbubble events, similar to the setup invoked by Dubois & Teyssier (2008). Our equilibrium setup allows us to investigate the reaction of the system to systematic changes of parameters like the halo pressure or the superbubble size. Our methods described in sect. 3 comprise the most important physics, however, some simplifications had to be made which require further discussion.
Firstly, our SN bubbles were triggered randomly, in accordance to the Kennicutt-Schmidt law for star formation (15), which is correlated to the column density , but neglects the volume density of the constituent cells within the column. One might argue that it would be better to let the SN probability increase with the density. This in turn will mean that more bubbles occur deeper within the disc and will thus have a harder time reaching the halo. The net effect would be an overall mitigation of the wind by an unknown factor. In consequence, energy would be converted even less efficiently.
We further had to implement a lower volume density threshold for cells to count as part of the disc and to amount to the surface density of their specific column. Note that the value used for our models, , is just a crude estimate for the lowest density regions found in the gas phase of the ISM, and thus allows for some variation. For instance, a lower threshold will open up a regime of rarefied cells surrounding the disc as is currently defined. This will have an effect on the distribution of the SN events, allowing for a bubble to blow out into the halo with less resistance. Though, the change in total will likely be of little effect regarding the wind strength - note that such rarefied cells will likely contain just around . This would definitely call for the modification of our probability function, which, if applied, would make an event in these cells extremely unlikely.
The star formation in our models is determined by a local Kennicutt-Schmidt law, and converges with increasing resolution. It is also sensitive to local events, such as material ejections and the bubble size. The latter is clearly a significant effect: In our resolution study (Figure 5), we find that the number of stars formed after agrees within 26 per cent. However, there is a convergence for resolutions finer than 36 pc. If we disregard the 65 pc resolution, the deviation already shrinks to nine per cent. Varying the bubble size yields a change in star formation of about 24 per cent. Larger bubble size leads to stronger star formation, yet very large bubbles seem to lead to such a strong outflow that the star formation gets weaker again. We believe that this feature of the model is realistic.
We have also investigated the dependence of the outflow rates on resolution. Apart from the 65 pc case (R65), which also shows a stronger deviation in the star formation rate, there is no strong or systematic deviation among all the runs with better resolution. One might argue that the smaller event sizes lead to less well-resolved superbubbles, and thus the higher outflow rates for larger superbubbles might be explained by such a resolution effect. The effect is however not present for the ST runs. Because the pressure in the disc is the same for ST and NC runs, the initial bubble sizes should be similar, which make resolution effects an unlikely explanation for the trends with bubble size.
Our simulations assume large-scale rotational symmetry. This means that the structure of the ISM is considerably restricted. Clumping in the azimuthal direction or spiral structure may increase the local star formation density, while leaving the average value constant. In a gas rich galaxy such as considered here local feedback effects should then also become important (e.g. Hopkins et al. 2012). In general, one should expect that superbubbles would break out of the dense gas more easily, as the effective surface of the star-forming regions is increased. This should reduce the critical bubble size for winds to occur and might increase the outflow rates. For these reasons our results should be regarded as qualitative, only.
5.2 Wind drivers
In our simulations, winds form from sequences of superbubbles: Individual superbubbles break out of the disc and remain in the lower halo for about yrs. We find hardly any bubbles to be affected by buoyancy. Either bubbles simply collapse after having overexpanded and fall back to the disc, or they are fed by subsequent superbubbles which eventually drives a larger scale outflow. For our standard run, this happens above a critical star formation density of 0.1 yr kpc. We show that the effect as such may happen even if we inject either only the thermal energy or only the kinetic energy component. The outflow is then however attenuated. There is a good reason why we should expect such a behaviour: A pressure supported bubble will simply expand into the direction of the strongest pressure decline, i.e. radially outward, once the halo is reached. Gas which just has its kinetic energy may not accelerate as efficiently by this pressure gradient. On the contrary the pressure force works the other way, because the bubble is under-pressured much faster. Moreover, we find the outflow energy to vary by a factor of 100 for different event sizes (Figure 17 and Table LABEL:tab-flux), while at the same time NC200 produces a strong wind, and NC100 a weak one, just strong enough to enter the thermally pressured halo. Hence, only superbubbles comprising at least around 100 SNe seem useful in order to generate a wind, slightly more if the halo is still relatively massive but cools inefficiently, as modelled in our NC runs.
Furthermore, in sect. 4.4 we have seen that the bubble size matters during the phase where the wind is launched and breaching through the inner halo regions. In reality, bubbles will not be all of the same size but rather occur in a wide range from single, isolated SNe to a few hundred per bubble. Here, we show that the larger superbubbles matter the most for galactic winds. However, for more realistic event size distributions the mass and energy flux in the resulting wind might converge earlier and exhibit fewer and smaller peaks. For the time being, we will leave this matter open for future investigation.
The bubble size has turned out to be relevant firstly for the initial shock wave, and secondly in the steady wind phase. Larger bubbles give a more powerful rise to the wind, and will keep their strength at higher, roughly constant levels for a long time. In LBGs which presumably blow winds continuously at a steady level, the bubble size could be important. Moreover, during a starburst phase, where the wind is often young and after which star formation will decrease rapidly, the bubble size might play a role.
The thermal halo pressure determines whether or not a stable wind phase develops in the first place. If the pressure is too high, the wind may very well be unable to proceed too far from the disc. How far it can go depends on the bubble size. Winds set up by small bubbles will stop early and in some cases collapse back onto the disc entirely, whereas winds resulting from large superbubbles have a good chance of escaping the halo no matter the halo pressure.
5.3 Comparison to observations
We model galaxies with masses, star formation rates, and star formation densities comparable to LBGs. We find that our superbubble mechanism starts to produce outflows from star formation densities around 0.1 yr kpc, which would predict outflows in essentially all LBGs, as observed. While we do not follow the ionisation state of the gas, it is clear that the gas around 1 cm in our simulations is highly turbulent, often escaping the disc (Figure 18). This is in good agreement with the generally observed high level of turbulence in LBGs, e.g. the metal line observations in Law et al. (2012). The superbubble mechanism produces naturally weak, radiative shocks in the halo, which we observe directly in the beginning of the simulations (Figure 19).
These shells would plausibly be observed as Ly absorption systems with column densities around cm to cm, and velocities around 200 km s, which compares favourably to observations (e.g. Verhamme et al. 2008). If anything, one would like to increase the halo density by a factor of a few, in order to cover the full range of observations.
At later times, these shells leave the computational domain. In many of our simulated galaxies, the shells actually come back and may start over. This happens when we moderate the feedback mechanism or superbubble size. Overall, the mechanism seems well able to explain wind shells in LBGs.
In observations of nearby galactic winds, one frequently finds energy efficiencies of order ten per cent (Veilleux et al. 2005). We find much less in our simulations. A similar discrepancy is seen in the mass outflow rates. It is well possible that these are different classes of objects. If the high wind efficiencies would also be confirmed for LBGs, it might point to some effects we might still be missing in our simulations. In particular, if the galactic halo is cleared out in some regions, simulations which more accurately treat the feedback mechanisms in the disc and which find much higher mass loss rates from the disc might become more applicable.
We have performed hydrodynamic simulations in spherical wedge geometry with polar cutouts with the grid-based 3D code NIRVANA, setting up a disc-halo system close to hydrodynamic equilibrium. Stars are formed in agreement to a local Kennicutt-Schmidt law, and superbubbles are injected according to a preset number of instantly occurring SNe. Thus, the structure of the interstellar medium is simplified to large-scale effects. Spiral-arm structure is not taken into account. In regions of the galactic discs, where a critical star formation density is exceeded, series of superbubbles expand into the galactic halo at slightly supersonic speed, lead by radiative shocks which may produce cold gas shells.
We parametrise the outflow strength via mass and energy loss rates. For our chosen galaxy setup and feedback implementation, the outflow strength depends on the details of the feedback implementation, the stellar content per superbubble and the halo pressure: We do not get significant outflows if we neglect the thermal energy part in the superbubble injection, and only retain the kinetic part. Outflows are stronger for bigger individual superbubbles, if the halo pressure is not too large. For large halo pressure, outflows are suppressed for superbubbles below a critical size. Our absolute outflow rates are uncertain due to the assumption of uniform stellar content per superbubble in a given simulation, and the simplifications in the modelling of the structure of the interstellar medium and the superbubble injection mechanism.
The simulated galaxies are in overall good agreement with LBGs regarding mass, star formation rates, star formation densities, gas kinematics and expanding shell characteristics. This suggests that the interaction between superbubble-driven winds and a heavy gaseous halo is a good candidate to explain the characteristics of LBGs.
The authors would like to thank the CAST group for the regular meetings including several useful suggestions on key aspects of this paper. We also would like to thank the Excellence Cluster Universe and the DFG (KR 2857/3-1) for the financial support, which made this project possible. Our cordial thanks further goes to Barbara Ercolano, Mark Westmoquette, Anne Verhamme and Daniel Schaerer for helpful discussions. We are also grateful to Ralf-Jürgen Dettmar for fruitful discussions on the topic during his visit to USM. Special thanks finally goes to Christian Tapken and Matthew Lehnert for helpful comments on our project during the DFG conference in Potsdam, September 2010. We also thank the anonymous referee for very useful comments that considerably improved this paper.
- Adelberger et al. (2003) Adelberger, K. L., Steidel, C. C., Shapley, A. E., & Pettini, M. 2003, ApJ, 584, 45
- de Avillez & Breitschwerdt (2004) de Avillez, M. A., & Breitschwerdt, D. 2004, A&A, 425, 899
- de Avillez & Breitschwerdt (2005) de Avillez, M. A., & Breitschwerdt, D. 2005, A&A, 436, 585
- Bagetakos et al. (2011) Bagetakos, I., Brinks, E., Walter, F., et al. 2011, AJ, 141, 23
- Barai et al. (2013) Barai, P., Viel, M., Borgani, S., et al. 2013, MNRAS, 430, 3213
- Bendo et al. (2010) Bendo, G. J., et al. 2010, MNRAS, 402, 1409
- Burkert (2006) Burkert, A. 2006, Comptes Rendus Physique, 7, 433
- Churazov et al. (2001) Churazov, E., Brüggen, M., Kaiser, C. R., Böhringer, H., & Forman, W. 2001, ApJ, 554, 261
- Cioffi et al. (1988) Cioffi, D. F., McKee, C. F., & Bertschinger, E. 1988, ApJ, 334, 252
- Cooper et al. (2008) Cooper, J. L., Bicknell, G. V., Sutherland, R. S., & Bland-Hawthorn, J. 2008, ApJ, 674, 157
- Creasey et al. (2011) Creasey, P., Theuns, T., Bower, R. G., & Lacey, C. G. 2011, MNRAS, 415, 3706
- Creasey et al. (2013) Creasey, P., Theuns, T., & Bower, R. G. 2013, MNRAS, 429, 1922
- Dalla Vecchia & Schaye (2012) Dalla Vecchia, C., & Schaye, J. 2012, MNRAS, 426, 140
- Dib et al. (2006) Dib, S., Bell, E., & Burkert, A. 2006, ApJ, 638, 797
- Dubois & Teyssier (2008) Dubois, Y., & Teyssier, R. 2008, A&A, 477, 79
- Erb et al. (2006) Erb, D. K., Steidel, C. C., Shapley, A. E., et al. 2006, ApJ, 647, 128
- Flynn et al. (1996) Flynn, C., Sommer-Larsen, J., & Christensen, P. R. 1996, MNRAS, 281, 1027
- Forman et al. (2007) Forman, W., Jones, C., Churazov, E., et al. 2007, ApJ, 665, 1057
- Franx et al. (1997) Franx, M., Illingworth, G. D., Kelson, D. D., van Dokkum, P. G., & Tran, K.-V. 1997, ApJL, 486, L75
- Frye et al. (2002) Frye, B., Broadhurst, T., & Benítez, N. 2002, ApJ, 568, 558
- Genzel et al. (2008) Genzel, R., et al. 2008, ApJ, 687, 59
- Gnedin (1998) Gnedin, N. Y. 1998, MNRAS, 294, 407
- Gourgouliatos & Lyutikov (2012) Gourgouliatos, K. N., & Lyutikov, M. 2012, MNRAS, 420, 505
- Greve & Sommer-Larsen (2008) Greve, T. R., & Sommer-Larsen, J. 2008, A&A, 480, 335
- Gupta et al. (2012) Gupta, A., Mathur, S., Krongold, Y., Nicastro, F., & Galeazzi, M. 2012, ApJL, 756, L8
- Heath et al. (2007) Heath, D., Krause, M., & Alexander, P. 2007, MNRAS, 374, 787
- Heckman et al. (1998) Heckman, T. M., Robert, C., Leitherer, C., Garnett, D. R., & van der Rydt, F. 1998, ApJ, 503, 646
- Heckman et al. (2000) Heckman, T. M., Lehnert, M. D., Strickland, D. K., & Armus, L. 2000, ApJS, 129, 493
- Hopkins et al. (2012) Hopkins, P. F., Quataert, E., & Murray, N. 2012, MNRAS, 421, 3522
- Kaiser et al. (2005) Kaiser, C. R., Pavlovski, G., Pope, E. C. D., & Fangohr, H. 2005, MNRAS, 359, 493
- Kennicutt (1998) Kennicutt, R. C., Jr. 1998, ApJ, 498, 541
- Kirkpatrick et al. (2011) Kirkpatrick, C. C., McNamara, B. R., & Cavagnolo, K. W. 2011, ApJ, 731, L23
- Krause (2005) Krause, M. 2005, A&A, 436, 845
- Krause et al. (2013) Krause, M., Fierlinger, K., Diehl, R., et al. 2013, A&A, 550, A49
- Kulas et al. (2012) Kulas, K. R., Shapley, A. E., Kollmeier, J. A., et al. 2012, ApJ, 745, 33
- Laursen et al. (2009) Laursen, P., Razoumov, A. O., & Sommer-Larsen, J. 2009, ApJ, 696, 853
- Law et al. (2012) Law, D. R., Steidel, C. C., Shapley, A. E., et al. 2012, ApJ, 759, 29
- Lowenthal et al. (1997) Lowenthal, J. D., et al. 1997, ApJ, 481, 673
- Mac Low et al. (1989) Mac Low, M.-M., McCray, R., & Norman, M. L. 1989, ApJ, 337, 141
- Mac Low et al. (1998) Mac Low, M.-M., Klessen, R. S., Burkert, A., & Smith, M. D. 1998, Physical Review Letters, 80, 2754
- Martin (2005) Martin, C. L. 2005, ApJ, 621, 227
- Miyamoto & Nagai (1975) Miyamoto, M., & Nagai, R. 1975, PASJ, 27, 533
- Moster et al. (2013) Moster, B. P., Naab, T., & White, S. D. M. 2013, MNRAS, 428, 3121
- Nath & Trentham (1997) Nath, B. B., & Trentham, N. 1997, MNRAS, 291, 505
- Navarro et al. (1996) Navarro, J. F., Frenk, C. S., & White, S. D. M. 1996, ApJ, 462, 563
- Ntormousi & Sommer-Larsen (2010) Ntormousi, E., & Sommer-Larsen, J. 2010, MNRAS, 409, 1049
- Oey & García-Segura (2004) Oey, M. S., & García-Segura, G. 2004, ApJ, 613, 302
- Oey (2009) Oey, M. S. 2009, American Institute of Physics Conference Series, 1156, 295
- Ohta et al. (2003) Ohta, Y., Kayo, I., & Taruya, A. 2003, ApJ, 589, 1
- Pettini et al. (2000) Pettini, M., Steidel, C. C., Adelberger, K. L., Dickinson, M., & Giavalisco, M. 2000, ApJ, 528, 96
- Pettini et al. (2001) Pettini, M., Shapley, A. E., Steidel, C. C., Cuby, J.-G., Dickinson, M., Moorwood, A. F. M., Adelberger, K. L., & Giavalisco, M. 2001, ApJ, 554, 981
- Pettini et al. (2002) Pettini, M., Rix, S. A., Steidel, C. C., Adelberger, K. L., Hunt, M. P., & Shapley, A. E. 2002, ApJ, 569, 742
- Powell et al. (2011) Powell, L. C., Slyz, A., & Devriendt, J. 2011, MNRAS, 414, 3671
- Putman et al. (2012) Putman, M. E., Peek, J. E. G., & Joung, M. R. 2012, ARA&A, 50, 491
- Rasera & Teyssier (2006) Rasera, Y., & Teyssier, R. 2006, A&A, 445, 1
- Razoumov & Sommer-Larsen (2007) Razoumov, A. O., & Sommer-Larsen, J. 2007, ApJ, 668, 674
- Richling (2003) Richling, S. 2003, Ap&SS, 284, 361
- Roediger et al. (2007) Roediger, E., Brüggen, M., Rebusco, P., Böhringer, H., & Churazov, E. 2007, MNRAS, 375, 15
- Rupke et al. (2002) Rupke, D. S., Veilleux, S., & Sanders, D. B. 2002, ApJ, 570, 588
- Rupke et al. (2005) Rupke, D. S., Veilleux, S., & Sanders, D. B. 2005, ApJS, 160, 115
- Scannapieco & Brüggen (2010) Scannapieco, E., & Brüggen, M. 2010, MNRAS, 405, 1634
- Schaerer et al. (2011) Schaerer, D., Hayes, M., Verhamme, A., & Teyssier, R. 2011, A&A, 531, A12
- Shapley et al. (2003) Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, ApJ, 588, 65
- Sharp & Bland-Hawthorn (2010) Sharp, R. G., & Bland-Hawthorn, J. 2010, ApJ, 711, 818
- Silk et al. (2012) Silk, J., Antonuccio-Delogu, V., Dubois, Y., et al. 2012, A&A, 545, L11
- Smail et al. (2003) Smail, I., Chapman, S. C., Ivison, R. J., Blain, A. W., Takata, T., Heckman, T. M., Dunlop, J. S., & Sekiguchi, K. 2003, MNRAS, 342, 1185
- Sommer-Larsen (2006) Sommer-Larsen, J. 2006, ApJL, 644, L1
- Sommer-Larsen & Toft (2010) Sommer-Larsen, J., & Toft, S. 2010, ApJ, 721, 1755
- Springel & Hernquist (2003) Springel, V., & Hernquist, L. 2003, MNRAS, 339, 289
- Steidel et al. (2010) Steidel, C. C., Erb, D. K., Shapley, A. E., et al. 2010, ApJ, 717, 289
- Stinson et al. (2006) Stinson, G., Seth, A., Katz, N., Wadsley, J., Governato, F., & Quinn, T. 2006, MNRAS, 373, 1074
- Stil et al. (2009) Stil, J., Wityk, N., Ouyed, R., & Taylor, A. R. 2009, ApJ, 701, 330
- Strickland & Stevens (2000) Strickland, D. K., & Stevens, I. R. 2000, MNRAS, 314, 511
- Sutherland & Dopita (1993) Sutherland, R. S., & Dopita, M. A. 1993, ApJS, 88, 253
- Swinbank et al. (2005) Swinbank, A. M., et al. 2005, MNRAS, 359, 401
- Tapken et al. (2007) Tapken, C., Appenzeller, I., Noll, S., et al. 2007, A&A, 467, 63
- Tegmark et al. (1993) Tegmark, M., Silk, J., & Evrard, A. 1993, ApJ, 417, 54
- Tremblay et al. (2012) Tremblay, G. R., O’Dea, C. P., Baum, S. A., et al. 2012, MNRAS, 424, 1026
- Veilleux et al. (2005) Veilleux, S., Cecil, G., & Bland-Hawthorn, J. 2005, ARA&A, 43, 769
- Verhamme et al. (2008) Verhamme, A., Schaerer, D., Atek, H., & Tapken, C. 2008, A&A, 491, 89
- Verhamme et al. (2012) Verhamme, A., Dubois, Y., Blaizot, J., et al. 2012, A&A, 546, A111
- Voss et al. (2009) Voss, R., Diehl, R., Hartmann, D. H., et al. 2009, AAP, 504, 531
- Weaver et al. (1977) Weaver, R., McCray, R., Castor, J., Shapiro, P., & Moore, R. 1977, ApJ, 218, 377
- Wünsch et al. (2011) Wünsch, R., Silich, S., Palouš, J., Tenorio-Tagle, G., & Muñoz-Tuñón, C. 2011, ApJ, 740, 75
- Zhao et al. (2009) Zhao, D. H., Jing, Y. P., Mo, H. J., Börner, G. 2009, ApJ, 707, 354
- Ziegler & Yorke (1997) Ziegler, U., & Yorke, H. W. 1997, Computer Physics Communications, 101, 54
Appendix A Control runs
We have performed control runs to check for the impact of our choice of the size of the azimuthal wedge and the polar cutout. Run STpol is identical to ST100 apart from the polar cutout having been reduced from to . Run STaz is again identical to run ST100, but the we simulate half the azimuthal angle, instead of just . We show a comparison of meridional density plots in Figure 20 at 10 Myr, when the wind first develops. From the plot, it is clear that the general structure of the disc ISM, and the extent of the superbubble expansion into the halo are very similar also to the corresponding plot for run ST100 (Figure 6). It is beyond the scope of this article to investigate the effect of prominent non-axisymmetric structures on outflow rates, like e.g. spiral arms.
The run of the outflow rates in run STpol are within a factor of a few similar to the ones of ST100. This is within the general scatter of the ST runs, and thus expected. We show the cumulative mass and energy fluxes over the polar angle in Figure 21, averaged over all the 201 snapshots of the simulation. In the standard runs, 0.8 per cent of the solid angle is cut out. In run STpol, only 0.2 per cent of the solid angle is still cut out. The extra 0.6 per cent of solid angle contribute 1 percent of the mass flux and 10 per cent of the energy flux. The differential energy flux declines towards the polar axis. This should be regarded as the typical uncertainty due to the polar cutout.