How runaway stars boost galactic outflows

Roughly ten per cent of OB stars are kicked out of their natal clusters before ending their life as supernovae. These so-called runaway stars can travel hundreds of parsecs into the low-density interstellar medium, where momentum and energy from stellar feedback is efficiently deposited. In this work we explore how this mechanism affects large scale properties of the galaxy, such as outflows. To do so we use a new model which treats OB stars and their associated feedback processes on a star-by-star basis. With this model we compare two hydrodynamical simulations of Milky Way-like galaxies, one where we include runaways, and one where we ignore them. Including runaway stars leads to twice as many supernovae explosions in regions with gas densities ranging from 1e-5 cm^-3 to 1e-3 cm^-3. This results in more efficient heating of the inter-arm regions, and drives strong galactic winds with mass loading factors boosted by up to one order of magnitude. These outflows produce a more massive and extended multi-phase circumgalactic medium, as well as a population of dense clouds in the halo. Conversely, since less energy and momentum is released in the dense star forming regions, the cold phase of the interstellar medium is less disturbed by feedback effects.


INTRODUCTION
Λ-cold dark matter has been highly successful in explaining and predicting a variety of observed properties, such as large scale structure, halo clustering and galaxy scaling relations (Eisenstein et al. 2005;Springel et al. 2005;Viel et al. 2008;Reid et al. 2010;Klypin et al. 2011;Komatsu et al. 2011;Somerville & Davé 2015). Nonetheless, this model has also encountered challenges related to the 'baryon cycle', i.e. how galaxies accrete and expel their gas. Galaxy formation is an inefficient process, with at most ∼ 1/3 of the cosmological baryon fraction being turned into stars in galaxies as massive as the Milky Way, and significantly less in dwarf galaxies (e.g. Behroozi et al. 2019). This notion has been notoriously difficult to predict by numerical simulations of galaxy formation, which historically have suffered from excessive gas cooling and loss of angular momentum, leading to simulated galaxies with little in common with observations (for a review on this topic, see Naab & Ostriker 2017).
Galactic outflows, driven by stellar feedback (Dekel & Silk 1986) and active galactic nuclei (Benson et al. 2003), are commonly believed to be instrumental for solving these problems (Somerville & Davé 2015). In the past decades, many studies have explored ways of numerically capturing the impact of stellar feedback processes on the interstellar medium (ISM), e.g. effects from supernovae (SNe) explosions, stellar winds, radiation from young stars (see, e.g., Katz 1992;Stinson et al. 2006;Agertz et al. 2013; E-mail: eric@astro.lu.se Keller et al. 2014;Simpson et al. 2015;Hopkins et al. 2018). This effort has made it possible to model galactic outflows in a cosmological context as an emergent property of clustered star formation, with simulations now matching a range of observed galaxy properties (e.g. Hopkins et al. 2014;Agertz & Kravtsov 2015, as well as properties of the turbulent ISM and giant molecular clouds (GMCs; Grisdale et al. 2017Grisdale et al. , 2018Grisdale et al. , 2019).
Adding to the success of stellar feedback models, simulations of entire galaxies has improved significantly in the recent years with parsec, or even sub-parsec spatial resolution and mass resolution reaching a few solar masses enabling models to capture the dense, cold ISM (see, e.g., Renaud et al. 2013;Rosdahl et al. 2017;Wheeler et al. 2018;Agertz et al. 2020). State-of-the-art simulations can now resolve most of the cooling radii of individual SN explosions which allows them to capture the build-up of momentum during the Sedov-Taylor phase which is crucial for robustly modelling SNe feedback (Kim & Ostriker 2015;Martizzi et al. 2015a). Furthermore, recent numerical studies have reached high enough resolution for individual stars to be modelled in galaxy scale simulations (Hu et al. 2017;Wheeler et al. 2018;Su et al. 2018;Emerick et al. 2019;Lahén et al. 2019). We stress however that despite the improvements in galaxy modelling in the past decade, (subgrid) stellar feedback still operates at the resolution level, with its coupling to the ISM not yet being fully understood (Ohlin et al. 2019).
While the actual process of star formation is not yet captured in galactic context, these simulations allow for star particles to represent individual stars, all sampled from an initial mass function (IMF) with subsequent feedback processes emerging on a star-by-star basis rather than from macro particles representing entire stellar populations. However, while galaxy simulations have started to resolve the stellar component star-by-star, they are still far from treating the stars collisionally. While a collisionless approximation is valid on galactic scales, star-star interactions drive the internal dynamics of dense star clusters. As a result, some fraction of OB stars are kicked out of their natal star clusters with velocities large enough to move them out of the dense star forming gas before exploding as SNe (e.g., Gies & Bolton 1986;Gies 1987;Stone 1991;Hoogerwerf et al. 2000;Schilbach & Röser 2008;Jilinski et al. 2010;Silva & Napiwotzki 2011;Maíz Apellániz et al. 2018). Simulation of clusters (e.g., Oh & Kroupa 2016) reveal that such velocity kicks originates from gravitational interactions (Poveda et al. 1967), including the disruption of binaries through SNe (Blaauw 1961). None of these effects can be taken into account self-consistently without collisional dynamics, which currently is computationally infeasible in galaxy scale simulations.
Massive runaway stars have been suspected to impact galaxy evolution through efficient injection of energy into low density gas (Naab & Ostriker 2017). Coupling energy from SNe to gas depends both on the phase of the gas (Katz 1992;Ceverino & Klypin 2009), as well as on the structure of the ISM (Kim & Ostriker 2015;Martizzi et al. 2015b;Walch & Naab 2015;Iffrig & Hennebelle 2015;Ohlin et al. 2019). Because runaway stars move away from their natal regions, SN energy from these objects ought to couple differently compared to the stars which remains in closer proximity to the dense star forming gas. Typically explosions in low density gas generate a hot bubble which expands and leaves the galaxy as an outflow (Korpi et al. 1999;de Avillez 2000;de Avillez & Breitschwerdt 2004;Joung & Mac Low 2006). A similar effect has also been studied by comparing feedback from SNe positioned at density peaks compared to random positioning (see e.g. Korpi et al. 1999;de Avillez 2000;de Avillez & Breitschwerdt 2004Joung & Mac Low 2006;Walch et al. 2015;Martizzi et al. 2016;Girichidis et al. 2016;Gatto et al. 2017). Although these experiments are idealised, targeting only a small patch of the galaxy, they demonstrate that feedback from randomly located SN tends to drive stronger outflows compared to explosions in density peaks. This potentially impact the way the host galaxy forms and evolves across cosmic time.
In this paper we present a new model for treating individual stars, focusing on the effect of runaways. The model has been applied to simulations of an isolated disc galaxy with properties similar to the Milky Way. We describe how we model individual stars as well as our treatment of runaway stars in Section 2. The numerical setup and results are described in Section 3 and 4 respectively, while we discuss their implications in Section 5. Finally we summarise in Section 6.

STAR FORMATION AND STELLAR FEEDBACK
In order to probe the effects of runaway stars, we implement a new recipe for treating massive stars as individual particles which evolves during the simulation. In this section we describe the method used for sampling stars from an IMF, how feedback via winds and SNe is treated and how it chemically enriches the ISM.

Sampling the initial mass function
Treating massive stars as individual particles requires a method to sample stellar masses. To do this efficiently we employ the method by Sormani et al. (2017). In a mass bin indexed by i, a number of stars n i is determined through realisations of a Poisson distribution given by where the mean of the distribution, λ i , is given for each bin by where for bin i, in which stars have an average mass of m i , there is a fraction of mass f i out of the total mass M, available for star formation. The total stellar mass generated in N bins is then which is approximately equal to the mass M, due to the stochastic sampling, but converges towards M for large numbers of stars.
There are two properties of this method which makes it ideal for our simulations. Firstly, the computational expense is determined by the number of adopted mass bins, since the method only samples one random number for each bin and not for every star. Secondly, the bin-sizes can be chosen arbitrarily, which allows for stars below a certain mass threshold to be grouped in a single bin, while stars above this mass are sampled at high mass resolution. A problem with this method is the non-zero probability of sampling a population of stars with mass larger than the available gas mass. To avoid this, we sample the IMF starting from the lowmass end until either reaching the most massive stars defined by the mass bins or until the available gas is depleted. In the first case we return the extra mass to the gas parcel where it was initially collected, while in the second case we remove the most massive stars such that the stellar mass formed does not exceed the initial gas mass. This results in an IMF which is slightly more bottom heavy than the one originally targeted, which in principle could lead to weaken stellar feedback. Nonetheless, this effect is negligible compared to the uncertainties in the models for stellar feedback.
Our model uses of the three part IMF from Kroupa (2001), where the number of stars of mass m is given by where A is a normalisation factor and C j are constants that ensures continuity in the intersections given by C 1 = 1, . The three different parts are split into mass regimes given by m 1 = 0.01 M , m 2 = 0.08 M , m 3 = 0.5 M and m 4 = M max , with slopes α 1 = 0.3, α 2 = 1.3 and α 3 = 2.3. The normalisation is determined for stars with masses between M min = 0.01 M and M max = 120 M .
In this work we sample stars from the aforementioned IMF at each star formation event. We divide stars into high mass stars (HMS) and low mass stars (LMS). The HMS are defined to have stellar masses larger than 8 M , sampled up to 40 M using 100 equisized bins (giving a mass resolution of ∼ 0.3 M ). We do not include stars more massive than 40 M in our model. Such stars are both exceedingly rare and their short lifetimes means that the distance traveled by very massive runaway stars is short. However, these are extremely luminous important sources of feedback, especially in regards to stellar winds. Our model for stellar winds (described in detail in Section 2.3) has a strong scaling with stellar mass, and including these most massive stars can result in too strong early feedback. Furthermore, we assume core-collapse SNe for all HMS, which is unrealistic for too massive stars. In the code, all HMS are treated individual particles. Since we are interested in the effect of feedback from massive (> 8 M ) runaway stars, the rest of the IMF, defining the LMS, are lumped together into a single macro particle to reduce the computational cost.

Runaway stars
In the simulations, particles are treated in a collisionless manner. To simulate runaway stars, we give velocity kicks to all particles representing individual stars (i.e. ≥ 8 M ) at their formation. The collisional effects leading to kicks operate on spatial scales many orders of magnitudes below the gravitational softening length of our simulations, warranting a 'sub-grid' approach.
The method assumes kicks distributed isotropically with velocities, v, sampled from a power-law distribution, where the β is a free parameter in our model. The value of β depends on several factors. Dynamical scattering events which generate kicks typically involve the interaction between a hard binary and a third star. The frequency of such encounters will therefore not only depend on the relaxation time of the cluster, but must also be sensitive to the binary fraction. Moreover, as mentioned earlier, the kicks can also originate from the disruption of binary system caused by the SNe of one of the two components. This corroborates that binary fraction is a important parameter of the distribution of velocities. Furthermore, these processes implies a delay between the formation of a star and the time when it gets a kick. Because of the complexity of this problem, estimates of β typically demands the use of N-body simulations (see, e.g., Eldridge et al. 2011;Perets & Šubr 2012;Oh & Kroupa 2016;Renzo et al. 2019), but estimates from observations also exists (e.g., Hoogerwerf et al. 2000Hoogerwerf et al. , 2001. Furthermore, Banerjee et al. (2012) found that the velocity distribution shows some dependence on the mass of the runaway stars, with more massive stars reaching higher velocities (in agreement with Oh & Kroupa 2016).
In this work, we choose a value of β = 1.8 and normalise the distribution for values between 3 km s −1 and 385 km s −1 , which we apply to all HMS particles at formation without adding any time delay. This leads to ∼ 14% stars with velocities > 30 km s −1 . This choice is motivated by model MS3UQ_SP in Oh & Kroupa (2016), corresponding to runaways from a clusters with a mass of 10 3.5 M and half mass radii of 0.3 pc. Early observations indicated values for the runaway fraction of 30% (Stone 1991), however this is large compared to more recent work. As noted by Maíz Apellániz et al. 2018 (see also Silva & Napiwotzki 2011), the results of Stone overestimates this fraction. Maíz Apellániz et al. found observational evidence that 10−12% of O stars and a few percent of B stars are so called runaway stars with significant peculiar velocities (> 30 km s −1 ) with respect to their natal environment (in agreement with models by Eldridge et al. 2011;Renzo et al. 2019). Furthermore, the sampled velocities applied to the stars will change through gravitational forces acting on the stars throughout their lifetime. This change depends on the local gravitational field and the entire galactic potential for stars with long enough life times. This sensitivity to environment, together with the difficulty of comparing the observed population of runaway stars to the un-evolved velocity distribution, make a simple universal model (e.g. Equation 5) uncertain. The reader should be cautious of this and note that the results we present could overestimate the effect of runaway stars.
As a first approximation, one can naively compute the distance stars travel before exploding as SNe by multiplying their mean velocity with the average main-sequence life time in the appropriate mass range. Using the velocity distribution given by Equation 5 we find that stars in the mass range 8 M to 40 M with Z = Z 1 travels roughly 350 pc before exploding as supernova. This simple model does not take into account deceleration from the gravitational field of the cluster from which it escapes, and is therefore a upper limit on the travel distance. Nonetheless it is of the same order of magnitude as more detailed studies (see e.g. Eldridge et al. 2011;Renzo et al. 2019). The distance is significant, being an order of magnitude larger than the average size of star forming molecular clouds in the Milky Way (e.g. Heyer et al. 2009), and of the same order of magnitude as the scale height of the cold gas disc.

Stellar feedback
HMS particles are considered as individually evolving stars and we consider mass loss, enrichment, momentum-and energy-injection from fast winds and core-collapse SNe. In massive stars (≥ 8 M ), radiation pressure is significant enough to push away their outer envelopes, giving rise to strong stellar winds for their entire mainsequence evolution (see, e.g., Willis & Garmany 1987;Cassinelli & Lamers 1987). For this work we employ a model with a mass-loss rate based on Dale & Bonnell (2008) given by where M b is the birth mass of the star. The scaling to the metallicity 2 was added to the mass-loss rate to account for the metallicity dependence of the photon coupling to the stellar envelope, which drives the wind. The metallicity exponent γ has been shown to range between 0.5 and 0.8, (Kudritzki et al. 1987;Vink et al. 2001;Mokiem et al. 2007;Vink 2011) and in this work we adopt γ = 0.5. The velocity of the fast winds typically ranges between 1000 km s −1 and 3000 km s −1 in the literature (see, e.g., Leitherer et al. 1992;Lamers & Cassinelli 1999). We use a value v w = 1000 km s −1 for our model, and inject the wind into the surrounding gas as a continuous source of momentum during the star's main sequence life time. HMS stars explode as core-collapse SNe at the end of their main sequence. The main sequence lifetime, t MS , for a star given its mass and metallicity Z (here expressed as the total mass fraction in elements heavier than He). To determine t MS our model uses a stellar age-mass-metallicity fit by Raiteri et al. (1996), who found where the coefficients are given by Note that this gives very similar main sequence lifetimes compared to the single-star evolution (SSE) formulae by Hurley et al. (2000).
When a star explodes we deposit 10 51 erg of energy into the gas at the location of the star particle. If the SNe cooling radius is not resolved, hence the momentum built up during the Sedov-Taylor phase is not captured self-consistently, we explicitly inject the momentum from this phase following the scheme by Kim & Ostriker (2015, see also Martizzi et al. 2015a). Our implementation is identical to that of  and Rhodin et al. (2019) and we refer the reader to these works for more details.
For LMS particles we consider mass loss, enrichment, momentum-and energy-injection for type Ia SNe and asymptotic giant branch (AGB) winds. Details of the scheme can be found in Agertz et al. (2013), with the adopted Kroupa IMF (Equation 4) truncated at at the LMS upper mass limit (8 M ).

NUMERICAL SETUP
We have implemented the method described above in the N-body + Adaptive Mesh Refinement (AMR) code RAMSES (Teyssier 2002). RAMSES treats stars and dark matter through collisionless dynamics using the particle-mesh method (Hockney & Eastwood 1981;Klypin & Shandarin 1983), with the gravitational potential calculated by solving the Poisson equation using the multi-grid method (Guillet & Teyssier 2011) for all refinement levels. The fluid dynamics of the gas is calculated using a second-order unsplit Godunov method with an ideal mono-atomic gas with adiabatic index γ = 5/3. The code accounts for metallicity dependent cooling by using the tabulated cooling functions of Sutherland & Dopita (1993) for gas temperatures of 10 4−8.5 K, and rates from Rosen & Bregman (1995) for cooling down to lower temperatures.
We have carried out two simulations of an isolated star forming galaxy: one where the high-mass single stars are initiated with the velocity kick from the velocity distribution given by Equation 5 using β = 1.4 (referred to as runaways), and a control simulation where high-mass single stars are not assigned any velocity kick (referred to as no runaways) as in traditional galaxy simulations.
The initial conditions are those of the isolated disc galaxy in the AGORA project (Kim et al. 2014(Kim et al. , 2016, set up to approximate a Milky Way-like galaxy following the methods described in Hernquist (1993) and Springel (2000). Briefly, the dark matter halo has a concentration parameter c = 10 and virial circular velocity v 200 = 150 km s −1 . This translates into a halo virial mass M 200 = 1.1 × 10 12 M within R 200 = 205 kpc. The total baryonic disc mass is M disc = 4.5 × 10 10 M with 20% in gas. The initial stellar and gaseous components follow exponential surface density profiles with scale lengths r d = 3.4 kpc and scale heights h = 0.1r d . The bulge-to-disc mass ratio is 0.125. The bulge mass profile is that of Hernquist (1990) with scale-length 0.1r d . The halo and stellar disc are represented by 10 6 particles each, and the bulge consists of 10 5 particles.
The galaxies are simulated in a box with sides 600 kpc with adaptive mesh refinement allowing for a minimum cell size of 9 pc. The refinement strategy uses a quasi-Lagrangian approach which ensures a roughly constant number of particles (∼ 8) in each cell which reduces discreteness effects (Romeo et al. 2008). Furthermore, the AMR scheme splits cells into 8 new cells, each with mass m refine = 4014 M , when their sum of stellar and gas mass exceeds 8 × m refine . Star formation uses a standard procedure with the star formation rate density in cells exceeding a density threshold, ρ , given by where SF is the star formation efficiency per free-fall time, ρ g is the cell gas density and t ff = 3π/32Gρ is the local gas free-fall time (see Agertz et al. 2013, for details). For this work we adopt ρ = 100 cm −3 and sample star particles with a mass resolution of 500 M , from which we remove mass for the HMS sampling at each star formation event. For the simulations considered here we use a constant SF = 5%. Observationally, the efficiency per free-fall time averages 1% in Milky Way giant molecular clouds (GMCs) (Krumholz & Tan 2007), albeit with a large spread (Murray 2011; Lee et al. 2016).

RESULTS
We reemphasise that including runaway stars provides a means to relocate massive stars away from their dense formation sites before they end their lives in SNe, and our goal here is to explore how this process affects the evolution of a Milky Way-mass spiral galaxy. The two galaxies are evolved for 250 Myr, a time during which they both are actively star forming. In the first 80 Myr the galaxies adjust to their initial conditions, and in the case of the runaways simulation there is a suppression in star formation after 120 Myr (approximately one orbital time at the disc scale radius). After these 120 Myr both simulations has stable star formation rate. In this work all results taken as averages over time exclude this early phase unless otherwise stated. Figure 1 shows projections along the line of sight for luminosity, gas density and temperature for the galaxies at time t = 250 Myr. Focusing first on the luminosity of the two galaxies, we find qualitative agreement between the two simulations, although we highlight that the runaways simulation has more stars in the inter-arm region as well as above and below the disc mid-plane. This indicates that the effect of runaway stars is subtle in the luminosity of the galaxy. However, the projected gas density and temperature largely differ between the two cases. Here runaway stars give rise to large (∼ kpc), hot (T > 10 6 K) bubbles in the inter-arm regions, as well as more prominent outflows (as seen on the edgeon views) compared to the no runaways model. In general, the gaseous halo is more structured with warm (T < 10 6 K) clouds found far away from the disc ( 5 kpc) in the runaways simulation.

Star formation rates and ISM properties
We next consider the star formation rate (SFR) of the two simulations, shown in the upper panel of Figure 2. Excluding the initial transient, both simulations have similar SFR (8−9 M yr −1 ) for the first 120 Myr, followed by a decrease to a SFR of 5 − 6 M yr −1 in the runaways simulation for the remaining time. This decrease is linked to the depletion of cold (T < 10 4 K) gas, which is shown in the lower panel of Figure 2. The simulation with runaway stars features less mass in the cold phase at all times, which ultimately leads to a decline in the SFR 3 .
During the period of high SFRs for both models ( Since the SFR is similar in both simulations at this epoch (thus having similar amounts of cold gas turning into stars), the runaways model must reduce the net gas mass more efficiently, either by removing it from this cold phase or inhibit warm gas from cooling more strongly compared to the no runaways model. To some extent both of these likely play a role. After 120 Myr, the cold gas mass loss rate changes to 7 − 10 M yr −1 in both simulations, with the runaways simulation now having less cold gas in total. For the runaways simulation this transition coincide with the decrease in SFR. This implies that even at low SFR, more cold gas is lost per unit of stellar mass formed in the runaways simulation.
To understand why the model with runaway stars predicts a more rapid removal of the cold ISM from the galaxy, we focus to the runaway mechanism, namely the effect of relocating SNe progenitors. Figure 3 shows the projected gas density for the cold ISM (T < 10 4 K). We selected outputs where both simulations included visible SN bubbles (t = 200 Myr). The propensity of feedback from runaway stars to create large bubbles in the inter-arm regions is evident. For the runaways simulation there are young massive stars present in the voids, whereas for the no runaways all massive stars overlap with the low mass star particles by design.
Initially, density contrasts are generated by stellar feedback and large scale gravitational instabilities, some of which evolves into bubbles. A crucial difference is that with the inclusion of runaway stars, the bubbles are resupplied with SN progenitors which keep injecting energy and momentum into them even though actual star formation is impossible in those regions. This repeated injection gives rise to kpc-scale bubbles in the runaways model, and as we will discuss later it leads to energy venting into the CGM. Throughout the runaways simulation many such bubbles form, albeit with seemingly random timing.
Quantitatively, we find a clear excess of SNe explosions in gas with 10 −5 cm −3 < ρ < 10 −3 cm −3 when including runaways, peaking at more than twice as many explosions at ρ ∼ 10 −4 cm −3   with respect to the no runaways case. This implies a direct injection of energy in low density regions in contrast with the typical injection sites of SNe. This mechanism is visible in Figure 3, and the resulting heating is visible in the face-on projection of the gas temperature in Figure 1. Finally, we note that the no runaways simulation features, not surprisingly, more SNe in gas with higher density compared to the runaways simulation.
Despite this seemingly violent impact on the cold ISM, the overall structure of the disc remains intact. To some degree, runaway stars alleviate the dense star forming gas from part of the feedback. This is demonstrated in Figure 4 were we show the structure of the cold gas in the disc. The left plot shows the vertical distance from the mid-plane of the disc within which half of the cold gas mass is located as function of radius. In the right plot we show the gas mass as function of vertical distance z. Although the introduc-tion of runaway stars cause large bubbles in the inter-arm gas, the overall structure of the cold ISM is very similar in both simulations. In fact, runaway stars produce a slightly thinner disc. Furthermore, the vertical velocity dispersion of the cold gas is σ z ∼ 20 km s −1 at radius R = 5 kpc for both simulations, with a roughly linear decrease to ∼ 5 km s −1 at 10 kpc after which it flattens out. This is in close agreement with observed HI kinematics in nearby spiral galaxies (Tamburro et al. 2009). Hence, while runaway stars drive strong turbulence and creates large bubbles it does not strongly alter the morphology of the cold disc. This allows for the co-existence of thin galactic stellar discs and strong galactic outflows, an otherwise challenging aspect of galaxy models (e.g. Agertz et al. 2011;Roškar et al. 2014).
In summary, runaway stars induce in more efficient coupling of feedback to diffuse gas. The more rapid depletion of the cold ISM found above is then a result of more vigorous galactic outflows which we quantify in the next section.

Galactic outflows
The top panel of Figure 5 shows the vertical outflow rate as a function of time for both models. The rates are computed by summing cell-by-cell contributions in a cylindrical region with a radius of 12 kpc centred on the galaxy for 4 kpc < |z| < 6 kpc, considering only gas flowing in the vertical direction away from the mid-plane. The galaxy with runaway stars features significantly stronger outflows at all times, and right after the onset of the galactic wind (t > 80 Myr) the mass loss rate is more than 10 times higher than without runaway stars. This strong outflows in the runaways simulation is likely linked to the reduction in cold gas mass, causing the reduced SFR. At 150 Myr this wind calms down, although it still removes mass at roughly twice the rate of the no runaways model. Furthermore, there are burst of outflows (the largest of which takes place at 200 Myr) linked to runaway stars generating a number of large (> kpc) bubbles (apparent in the middle panel in Figure 1). In the bottom panel of Figure 5 we present the evolution of the  Figure 6. Top: Mean vertical velocity of the gas as function of vertical height, with warm gas (2 × 10 4 < T < 10 6 K) as solid lines and hot gas (10 6 < T < 10 8 K) as dotted lines. Bottom: Ratio of outflow rates (green) and mass loading factor (orange) between the two simulations (runaways divided by no runaways) as function of vertical distance. To compute the outflow rates we used bin heights of 1 kpc for |z | < 10 kpc, 2.5 kpc for 10 kpc < |z | < 30 kpc and 5 kpc for |z | > 30 kpc in order to account for the decreasing resolution with increased vertical height. The lines shows the mean of all simulation outputs ignoring the first 80 Myr. mass loading factor, defined as Since η is very sensitive to the height where it is measured, we show it for three different |z| corresponding to roughly 0.025, 0.25 and 0.5 times the virial radius of the galaxy. For a given model, η varies up to three orders of magnitude depending on where it is measured. However, the ratio between the simulations is mostly independent of where we measure it. We will return to this notion later as well in this section, as compare our models to observations in Section 5.3. Focusing on the innermost measurement, we find that at 100 Myr, η reaches its maximum at ≈ 0.5 in the no runaways simulation. With the stronger outflows measured in the runaway model, combined with the reduced SFR, the mass loading factor is significantly higher throughout the entire duration of the simulation, with periods of η reaching values of ∼ 5 − 10. As we discuss further in Section 5.1, such an efficient wind driving implies a markedly different long term evolution for the galaxy.
Higher outflow rates do not necessarily imply faster galactic winds. In fact, in the runaways simulation we find that warm outflows (2 × 10 4 < T < 10 6 K) carry more mass but travel at a lower velocity (see solid line in bottom panel of Figure 6) than in the no runaways simulation. If the galactic winds were purely momentum-driven, the outflow mass would depend on the wind velocity as m ∝ 1/v. We measure the velocity to be less than twice as high in the runaways which then would imply an outflow mass no more than twice as high. However, the strong outflows in runaways yield ∼ 3 times more gas mass in the halo (as shown in Section 4.3). Therefore, this feedback is not momentum-driven. Conversely, in the case of purely energy-driven winds, one would expect m ∝ 1/v 2 , i.e. at most four times more in the runaways, which is compatible with our mass measurements.
In both models, v z increases with vertical distance, indicating that outflowing gas is accelerated by the continuously shock heated hot phase. Moreover, in contrast to no runaways, the runaways simulation features epochs of vigorous winds, with bursts of fast moving, hot (T > 10 6 K) gas originating from the large bubbles described earlier.
In the top panel of Figure 6 we compare the mean outflow properties between the two models as function of vertical distance after 120 Myr. Regardless of vertical distance from the disc, the model accounting for runaway stars always features an increased mass loss rate as well as mass loading factor. The impact of runaway stars on winds is hence not restricted to a narrowly defined region, but rather acts like a outflow strength 'multiplier' all the way up to ≈ 100 kpc. If generic, we note that this effect could allow for a simple treatment of the effect of runaway stars in semianalytical models of galaxy formation (see e.g. Somerville & Davé 2015) or cosmological simulations with phenomenological treatments of galactic winds (e.g. Vogelsberger et al. 2013).

Structure of gaseous halo
As already shown, runaway stars drive a strong galactic winds, launching large quantities of gas out in the halo. Despite the absence of cosmic inflows in our simulations and low resolution (∼ 100 pc) in the halo, the circumgalactic medium (CGM) produced by the outflows does play a crucial role in the long term evolution of the galaxy, e.g. by determining the recycling of gas used for star formation. Unless otherwise stated, we focus now on gas within 40 kpc of the galaxy centre from which we remove the disc (±2 kpc vertically and 15 kpc radially). Figure 7 shows the projected density shown for three different temperature ranges. Runaway stars produce a CGM which is more extended with the presence of long gas filaments and small clouds. The clouds rise from the disc and dissipate on timescales of tens of Myrs, in agreement with models of lifetimes of high velocity clouds around the Milky Way (Heitsch & Putman 2009). In our runaways simulation the motion of the clouds varies with some clouds being accreted back onto the disc, some buoyantly floating above the plane and some lifted outwards until dissipation. The filaments are the result of ram pressure exerted on these clouds by the hot medium.
At t = 250 Myr the gas mass in the considered halo region is 9.0 × 10 8 M and 2.9 × 10 8 M in the runaways and no runaways models respectively. The excess mass in the runaways model is in very hot (> 10 6 K), low density (< 10 −2 cm −3 ) gas, shown in the contours of Figure 8 (as well as the histograms). The background colour of this plot gives the mass ratio between runaways and no runaways in different phases, where we see that the runaways model gives over one order of magnitude more gas mass in this phase, demonstrating the more efficient coupling of SNe energy to the halo.
In the density histogram of Figure 8 we also see more mass for ρ ≈ 1 cm −3 in the runaways model. This is gas with temperature ∼ 10 4 K which sits in clouds and the filamentary structures. Finally, we note that in the no runaways simulation there is more diffuse (ρ ∼ 10 −5 cm −3 ) gas in the range 10 5 < T < 10 7 K than in the runaways simulation. runaways (T < 2 × 10 4 K) (2 × 10 4 < T < 10 6 K) (10 6 < T < 10 8 K) no runaways In all three phases the outflows caused by the inclusion runaway stars provide a significantly larger vertical scale. When comparing the time evolution of the two simulations in this view, we find repeated large bubbles (most clearly seen here in the hot gas) in the runaways simulation. These are absent in the no runaways simulation. The snapshot shown here is the last output (250 Myr).

Implications for long term galactic evolution
We have demonstrated that simulations incorporating massive runaway stars feature galactic outflows with higher mass loading factors compared to models without. This leads to a more massive and structured galactic halo, and an overall reduction of the mass in the cold ISM. The star forming ISM is not only regulated more efficiently by galactic outflows, some fraction of the runaway stars travel far enough to deposit energy into the inter-arm regions and the inner galactic halo. This increases the thermal energy in diffuse gas, which in turn reduces the gas accretion rate of the galaxy and thus acts as a preventive feedback. Such a process is energetically more efficient than mechanically ejecting the same amount of gas from the ISM. This means that stellar feedback with runaway stars does Ratio (runaway/no runaway) Figure 8. Comparison between the runaways (red) and no runaways (blue), as seen in the phase of the halo gas at the end of the simulation. The lines shows the mass distribution in gas phase, where contours are given for 1%, 10%, 40% and 80% and histograms shows the sum of all values along a given axis value. The background map shows the ratio between the models, coloured such that red indicated excess in runaways and blue indicated excess in no runaways.
not only lead to increased mass loading factors in galactic winds, but that it also includes a higher energy deposition into halo gas (Li & Bryan 2019).
Although our simulations have only been carried out for 250 Myr due to the high computational cost associated with treatment of individual stars, the above mechanisms imply a more efficient regulation of star formation over Gyr timescales, and hence an overall impact on galaxy evolution over a Hubble time. Applying our model in fully cosmological simulations, even for a limited amount of time, would be of great interest for shedding light on how runaway stars affect earlier phases of galaxy formation (z > 1) when gas accretion rates are significantly higher than at the current epoch and efficient stellar feedback is know to be important (e.g. . Runway star physics may therefore aid in explaining the inefficiency of galaxy formation from e.g. abundance matching (Behroozi et al. 2019), and will impact the precise galaxy mass range over which AGN feedback dominates stellar feedback (Benson et al. 2003).

Comparison to other work
The literature on galaxy scale simulations including runaway stars is limited due to the high computational cost of treating individual stars. Previous studies have therefore used alternative methods to overcome this challenge. Simulations of vertically stratified patches of discs have been used to compare clustered and random placement of SNe (see e.g. Korpi et al. 1999;de Avillez 2000;de Avillez & Breitschwerdt 2004 Early works studied SNe placement based on density thresholds (Korpi et al. 1999) and showed that some SN explosions must occur in isolation in order to produce a realistic multi-phase ISM (see, e.g., de Avillez 2000;de Avillez & Breitschwerdt 2004;Joung & Mac Low 2006;Walch et al. 2015;Li et al. 2017). In agreement with our results, these authors have found an increased effect from stellar feedback as a result of randomly placed SNe. In fact, models with stellar feedback injected over a broad range of densities leads to a very turbulent ISM and an overall reduction of star formation. Additionally, Kim & Ostriker (2018) found an increased mass loading factor as a result of including runaway stars, especially for hot gas, although this depends quantitatively on the vertical structure of the disc (Li et al. 2017).
The volume filling factor of gas with different temperatures has, not surprisingly, been shown to be affected by randomly placed SNe (see, e.g., de Avillez & Breitschwerdt 2004;Walch et al. 2015). This is also the case in our work, with more of the volume filled by hot gas (T > 3 × 10 5 K) in the disc when including runaway stars (71% vs. 54%). We find that the opposite is true for the halo, where the dense clouds expelled by the strong winds somewhat reduce the volume filling factor in hot gas.
For our model we account for the entire galactic disc, although this imposes constraints on the spatial resolution we can afford. As was pointed out by Martizzi et al. (2016), not accounting for the global geometry of the disc, as is the case for ISM patches, leads to unreliable wind properties due to incorrect treatment of the gravitational potential and the disc rotation.
To our knowledge, no other simulation work of entire disc galaxies have included treatment of individual runaway stars, al-though alternative methods have been studied. A model by Ceverino & Klypin (2009), also used in simulations by Ceverino et al. (2014) and Zolotov et al. (2015), includes the effect of runaway stars by applying velocity kicks to 10 − 30% of all stellar particles (with masses > 10 4 M ) representing entire stellar populations rather than individual massive stars. Ceverino et al. (2014) showed that their model depend strongly on resolution with its effect disappearing at low-resolution (60 pc) compared to high-resolution (14 pc) simulations. This is not unexpected since the mean travel distance of OB runaway stars is of order 100 pc (see Section 2.2) which sets a characteristic scale that needs to be resolved for these stars to differ from massive stars remaining in their natal regions.
In order to alleviate the computational expense required by tracing the trajectory of individual stars, an intermediate step could be to distribute SNe explosions around their natal stellar population. Tress et al. (2019) used this approach by injecting SNe energy at locations randomly sampled from a 5 pc Gaussian distribution. Although this captures some aspects like the spatial extent of star clusters, this is not representative of runaway stars which can travel significantly further and yield more complex distance distributions.

Comparisons to observations
Observational evidence for strong galactic winds is today ubiquitous e.g. in starburst galaxies (Veilleux et al. 2005;Rupke 2018). In this section we discuss how our simulations compare with the current observational data. We note that this comparison should be treated with caution since our simulations are idealised, isolated disc galaxies without their cosmological environment.
Mass loading factors in external galaxies have been observed to depend on galaxy mass although there is a considerable spread with mass loading factors ranging between η ∼ 0.01 − 10, even at a fixed stellar mass (Martin 1999;Bouché et al. 2012;Chisholm et al. 2017;Schroetter et al. 2019). Mass loading factors are difficult to extract from observations due the incomplete census of various gas phases as well as their sensitivity to where in the galaxy one measures them. Since observations typically use quasars absorption spectra to study winds in external galaxies the random positioning of the quasar with respect to star forming galaxies makes this a challenging endeavour. In Figure 9 we show a composite of observationally derived mass loading factors as function of stellar mass including a fit from Chisholm et al. (2017) together with our results. The data indicates a negative trend with stellar mass, as well as a large scatter. We show the final values of our results at three different heights above the disc plane. These values are in broad agreement with the entire spread in the observations, with the exception of the outermost measurements of the no runaways simulation. This implies that the spread in our simulations is due to either: 1) height above the disc plane; 2) the different stellar physics we consider, i.e., runaway stars. In the absence of information on the height of measurements (due to uncertain inclination of the disc) it is not possible to disentangle these two possibilities.
Other simulation work have modelled self-consistent driving of galactic outflows from stellar feedback in massive galaxies. Notably, Muratov et al. (2015) analysed the FIRE simulation suite and found that for M 10 10 M , mass loading factors were consistent with η ∼ 0 when measured at 0.25R vir ≈ 50 kpc (see their figure 6). This is incompatible with the observations compiled in Figure 9, and illustrates that the role of stellar feedback driven winds in galaxies of this mass is not yet understood. The results of Muratov et al. may indicate a lack of additional physics such as . Mass loading factor as function of stellar mass comparing the results of our simulations to values observed. Our simulations, with runaways in red and no runaways in blue, are shown measured at three different vertical heights (see Figure 5 for details). Orange markers show data for galaxies at redshift ∼ 0.1 from Bouché et al. (2012). Black diamonds shows data from nearby star-forming galaxies from Chisholm et al. (2017) with mass loading estimates computed from maximum outflow rates divided by star formation rates derived from UV spectra. The scaling relation derived in their work is shown as the dashed black line. The green crosses shows measurements from galaxies at redshift ∼ 0.1 by Schroetter et al. (2019). Although η depends strongly on the where outflows are measured the ratio between the two simulations remain the same. runaway stars as our no runaways simulation also produce low, albeit non-zero, mass loading factors.
Another striking effect of the runaway stars is the more structured and extended gas halo. The CGM is indeed observed to be highly structured (e.g. Werk et al. 2014, and references therein), with stellar feedback leaving a specific imprint on the properties of the circumgalactic gaseous haloes (e.g. Liang et al. 2016). To date, no consensus exists on what is driving the fine scale structure of the CGM, partly due to the fact that galaxy scale simulations are far from resolving the (likely) parsec scale clouds that make up the halo on small scales (see discussion in Section 5.5). We leave a study of the impact of run-away stars on halo absorber column density profiles and covering fractions for a future study.
In Section 4.3 we noted that more halo gas clouds existed across a wide range of temperatures in the runaways simulation. They arise from gas lifted from the disc as well as halo gas compression from self-gravity followed by efficient cooling (indeed, cooling is most efficient at T ∼ few × 10 5 K, with feedback driven perturbations leading to thermal instabilities Joung et al. (2012) and cloud formation). In reality, these clouds may be destroyed via hydrodynamical instabilities as they fall back onto the disc (Heitsch & Putman 2009), a process that is likely not fully captured in our models. Regardless, these clouds could be interpreted as high velocity clouds (HVC), which are observed in the Milky Way halo (Blitz et al. 1999;Wakker 2004;Moss et al. 2013Moss et al. , 2017. The origin of the HVCs is still an open question (Lockman et al. 2019), and although the most distant ones are likely from cosmic origin (see e.g., Richter 2012), there are plenty of clouds observed within few tens of kpc (see Wakker et al. 2008), compatible with our runaways model. It is likely that at least some fraction of these clouds originate from outflows (Quilis & Moore 2001), making the more effi-  Figure 10. Resolved Kennicutt-Schmidt relation, Σ SFR versus Σ g , comparing the runaway model (left) to the model ignoring runaway stars (right). Each point is a measurement from a 1 kpc square in the 24 kpc wide face-on view of the galaxy using data from our last output. To compute the star formation rate we use an age bin of 2 Myr. The colours show the mean velocity of all these stars, calculated after removing the velocity component coming from the rotational velocity. The grey dashed line shows the empirical relation fitted by Daddi et al. (2010) for disc galaxies. The runaways simulation includes a group of points with star formation rates ∼ 10 −5 M yr −1 kpc −2 which corresponds to runaway stars, as is indicated by the high velocity of the stars in these points. cient coupling of feedback energy to the galactic halo allowed by runaway stars an interesting addition to the formation scenario of HVCs.
Finally, we briefly highlight a curious effect found when including runaway stars in our simulations, namely their effect on the resolved Kennicutt-Schmidt (KS) relation. This relation is shown for the final output of both the runaways and no runaways simulation in Figure 10, together with the relation from Daddi et al. (2010) (see also Kennicutt 1998). With runaways we find regions of the disc which reside significantly below the canonical value, reaching Σ SFR ∼ 10 −6 M yr −1 kpc −2 . Since our model has a finite mass resolution for star formation events (∼ 500 M ), a lower limit to star formation surface densities of Σ SFR ≈ 10 −4 M yr −1 kpc −2 , for the accounted time (2 Myr) and area (1 kpc 2 ), is expected. This is indeed the case for no runaways, since stars does not significantly change in mass over 2 Myr. In contrast, in the runaways model this is enough for a fraction of runaway stars to venture away from their natal region and masquerade as inefficient star formation events. This is highlighted by the deviation in mean peculiar velocity of the stars shown by colour in Figure 10. Runaway star may therefore provide an explanation for the extremely inefficient star formation observed in local spiral and dwarf galaxies (Bigiel et al. 2010) which we will focus on in forthcoming work (Andersson et al., in prep.) 5.4 The impact of runaway stars in dwarf galaxies Dwarf galaxies have shallow potential wells and are expected to be strongly affected by stellar feedback, with higher mass loading factors than the cases studied here (see Figure 9) resulting in low galaxy formation efficiencies (Behroozi et al. 2019).
To give us an indication of how runaway stars impact these low mass galaxies we set up two additional pairs of simulations.
The first consists of a galaxy formed from collapsing gas in an analytical description of a dark matter halo. Its potential is set up according to a Navarro-Frenk-White (NFW) profile (Navarro et al. 1997) with R 200 = 50 kpc, where the rotational velocity is 35 km s −1 . This setup is identical to that of Teyssier et al. (2013), except they used a "live" halo described by particles. We adopt the same star formation and feedback recipes as is described in Section 2 and compared the evolution of one with and one without our runaway star model. After a settling phase the star formation and outflow rates were almost identical regardless of whether runaway stars are included or not. In addition to the isolated dwarf galaxy we simulated the formation of a ultra faint dwarf (UFD) galaxy (M 200 = 10 9 M at z = 0) in cosmological setting. The setup used is the same as the fiducial simulation in Agertz et al. (2020), except that we include our star formation and feedback models. Akin to the isolated dwarf galaxy we find a surprisingly small effect when including runaway stars, despite the mean travel distance for SNe progenitors being on the order of the half mass radius (∼ 300 pc) of the simulated UFD.
We attribute the lack of impact of runaways on dwarf galaxies, at least in terms of outflow properties, to the high porosity of the ISM which arises due to the shallow potential. This implies that energy from any SN is likely to efficiently couple to diffuse gas within or outside galaxy, regardless of the inclusion of runaways. Nonetheless, since runaway stars can leave the galaxy entirely one may expect a different enrichment history for the CGM. How this affects the chemistry of the dwarf galaxy and its environment in the long term deserves a dedicated study.

Limitations of the model
Throughout this work we have outlined limitations to our model. In this section we summarise these points, and discuss possible im-provements for the future. We highlight that this work is an exploratory study with the aim of giving a first approximation as to how runaway stars affect the physical processes in full galaxy simulations.
Treating self-consistently the star-star interactions responsible for the velocity kicks is still not possible in galactic context. This is due to the extreme computational cost a star-by-star description would represent, and also because of the need of a direct collisional treatment of gravity, to resolve close encounters and binary evolution, which is not possible with galaxy simulation codes. Therefore, we must rely on sub-grid models to assign the kick velocities, and ignore possible variation with the environment. Nonetheless, future work should focus on how the results here depend on the velocity distribution, which in turn can be studied outside the galacic context (see e.g. Eldridge et al. 2011;Perets & Šubr 2012;Moyano Loyola & Hurley 2013;Oh & Kroupa 2016;Renzo et al. 2019) Furthermore, the results shown here assumes a runaway fraction of 14% which is applied to all massive stars. Although this agrees with runaway fraction for O stars in a 10 3.5 M cluster (Oh & Kroupa 2016), it is not universal and observational data indicates that this fraction can be significantly lower with variations with stellar type. Maíz Apellániz et al. (2018) estimated the fraction at 10 − 12% for the O stars and ∼ 6% for B stars, which is also in agreement with those from Eldridge et al. (2011) when correcting for completeness. The results shown here could therefore overestimate the impact of runaway stars. Future work should focus on how the results depend on changing there runaway parameters.
The resolution of our simulations is limited to 9 pc. This is allows us to resolve the typical size of of the massive star forming gas clouds with approximately ten resolution elements, and this is essentially the scale runaway stars need to travel to reach low densities. Since this contrast in density is key to how the runaway stars affect the galactic evolution, an increase in resolution which allows for higher density contrast may also increasing the effect of runaway stars -in essence, less travel distance is required for more of the stars to reach low density gas to explode in. A suite of galaxy simulations with incrementally higher spatial resolution would be required to understand whether convergence can be reached.
Finally we note that due to the adopted adaptive refinement scheme the cooling length of the halo remains unresolved and we can only demonstrate systematic effects of the strong galactic wind caused by runaway stars.

CONCLUSIONS
Using hydrodynamical simulations of Milky Way-like galaxies, we have investigated how their evolution is affected by the inclusion of runaway stars, a mechanism not commonly accounted for in galaxy simulations. Our model initialises massive stars as individual particles with velocity kicks randomly sampled from a typical distribution for natal star clusters. For our star-by-star treatment, we have implemented a new model for stellar feedback accounting for fast stellar winds and core-collapse SN. We compare the impact of this implementation to one with similar feedback physics but without any velocity kicks, thus ignoring runaway stars. In summary we find the following differences when including runaway stars: (i) A significant fraction of SNe explode in low density gas (ρ < 10 −3 cm −3 ), with more than a doubling of SN explosions at densities 10 −4 cm −3 . By placing SNe progenitors in hot bubbles, where star formation is otherwise inhibited, the runaway mechanism significantly heats the inter-arm regions and drives strong galactic winds. This yields mass loading factors η ∼ 5 when measured at the disc-halo interface -an order of magnitude increase compared to models neglecting runaways. Both our models are compatible with observational values of η, which have large scatter due to uncertainty in where in the galaxy it is measured. Therefore, measurements of η with accurate estimates of location are required to fully understand the importance of accounting for runaway stars.
(ii) In the halo we find three times more gas mass, primarily in the hot diffuse phase, as well as a population of dense clouds which dissipate on timescales of tens of Myrs. This is mainly due to the large amount of gas lifted out by outflows. Furthermore, runaway stars are able to travel far enough to directly deposit SN energy into the halo, efficiently heating the CGM and preventing re-accretion of gas onto the disc.
(iii) The cold ISM (T < 10 4 K) is less disturbed by stellar feedback since runaways can leave the star formation regions. This implies that, although runaway stars lead to increased feedback efficiency, there is little change to the overall structure of cold ISM. Therefore, the effect of runaway stars cannot be simply modelled by increasing stellar feedback, but an actual relocation of the feedback sources must be implemented.
One limitation of our model is the assumption of a universal velocity distribution for stellar birth environments. Contrary to this, properties such as cluster density, binary fraction and IMF variations ought to play a role. Models accounting for stellar-scale properties (e.g. binary stars) in cosmological models exist (see e.g., Rosdahl et al. 2018), although without treatment for runaway stars, which requires tracing individual trajectories. A complete understanding of how runaways affects galaxy evolution might require a combination of these kind of models. However, treating even just a fraction of the stars as individual particles is very numerically expensive, therefore studying how our results depend on this velocity distribution is left to future work. Another limitation is that our model for stellar feedback ignores radiative processes. The strong luminosity of OB stars heats the surrounding medium before they explode as SN, thus affects how the energy is transferred to the gas. Adding radiative transfer to our already computationally heavy method is however beyond the scope of disc galaxy simulations. This would however be feasible for less massive galaxies, however, we do not see any significant effect of runaway stars for these galaxies. We attribute this to the high porosity of small galaxies implying that SNe location plays a minor role. As a final remark, the effect of runaway stars is clear in spiral galaxies for which the strong outflows lead to a more efficient regulation of star formation and in the long term this might impact the high end of the luminosity function.