The Structure and Composition of Multiphase Galactic Winds in a Large Magellanic Cloud Mass Simulated Galaxy

We present the first results from a high-resolution simulation with a focus on galactic wind driving for an isolated galaxy with a halo mass of ∼1011 M ⊙ (similar to the Large Magellanic Cloud) and a total gas mass of ∼6 × 108 M ⊙, resulting in ∼108 gas cells at ∼4 M ⊙ mass resolution. We adopt a resolved stellar feedback model with nonequilibrium cooling and heating, including photoelectric heating and photoionizing radiation, as well as supernovae, coupled to the second-order meshless finite-mass method for hydrodynamics. These features make this the largest resolved interstellar medium (ISM) galaxy model run to date. We find mean star formation rates around 0.05 M ⊙ yr−1 and evaluate typical time-averaged loading factors for mass (η M ∼ 1.0, in good agreement with recent observations) and energy (η E ∼ 0.01). The bulk of the mass of the wind is transported by the warm (T < 5 × 105 K) phase, while there is a similar amount of energy transported in the warm and the hot phases (T > 5 × 105 K). We find an average opening angle of 30° for the wind, decreasing with higher altitude above the midplane. The wind mass loading is decreasing (flat) for the warm (hot) phase as a function of the star formation surface rate density ΣSFR, while the energy loading shows inverted trends with ΣSFR, decreasing for the warm wind and increasing for the hot wind, although with very shallow slopes. These scalings are in good agreement with previous simulations of resolved wind driving in the multiphase ISM.


INTRODUCTION
Galactic winds are common across galaxy masses and sizes, and are generally believed to play a crucial role in regulating star formation and the galactic baryon budget, as well as the growth of supermassive black holes (SMBH).The role of galactic winds is twofold as they actively eject gas from the ISM, transporting mass, energy and metals into the circumgalactic medium (CGM) where they can then prevent the hot X-ray emitting gas in the CGM from falling victim to radiative cooling processes that would lead to over-cooling of the hot CGM gas.Thus, it has become apparent that galactic winds are important contributors to the assembly and formausteinwandel@flatironinstitute.orgtion of galaxies and represent one of the key aspects that needs to be understood to develop a detailed theory of galaxy formation and evolution.
The signatures of galactic winds are observed in the local universe (e.g.Lehnert & Heckman 1996;Martin 1999;McQuinn et al. 2019;Xu et al. 2022) as well as at high redshift (Pettini et al. 2001;Shapley et al. 2003;Förster Schreiber et al. 2019).Some observational work suggests that there is a rather weak correlation between outflow velocity and star formation rate (SFR) (Chisholm et al. 2015), while other groups suggest that there is a stronger correlation between outflow velocity and SFR (see e.g.Martin 2005;Heckman et al. 2015).
While the trend in the literature on the correlation between outflow velocity and galaxy mass is still somewhat under debate one can find more consensus in the literature on the mass loading factor η m = Ṁout /SFR, which tends to decrease with increasing galaxy mass which is not only established by observations at low redshift (e.g.Newman et al. 2012;Chisholm et al. 2017), but also appears to be robust against trends found at higher redshift (e.g.Förster Schreiber et al. 2014Schreiber et al. , 2019)).Similarly, one can find evidence that η m is positively correlated with star formation rate surface density (Σ sfr ), which is again a trend that manifests at lower redshift (e.g.McQuinn et al. 2019) and in high redshift systems (e.g.Davies et al. 2019).
However, despite the important findings on scaling relations of wind proprieties, galactic winds are much more complex than simple mass transport at a given outflow velocity, due to the observed multiphase nature of the winds.This fact makes observations of galactic winds extremely complicated since the cold, warm, and hot gas phases each require different observational techniques and/or instruments.The cold gas can be traced by neutral hydrogen (e.g.Walter et al. 2002;Martini et al. 2018) or via CO (e.g.Bolatto et al. 2013;Cicone et al. 2014;Leroy et al. 2015) and metals such as NaD, OI or MgI (e.g.Heckman et al. 2000;Rupke et al. 2005;Chen et al. 2010), and is characterized by temperatures between 10 and a few 100 K with rather low outflow velocities.The warm phase of the outflow, with temperatures around 10 4 K, can be studied by utilizing either Hα (e.g.McKeith et al. 1995;Westmoquette et al. 2009) or metal lines of photo-ionized species such as MgII, OIII and SiIV (e.g.Martin & Bouché 2009;Rubin et al. 2011;Nielsen et al. 2015;Chisholm et al. 2017).Even hotter winds (10 5 − 10 6 K) can be traced by the higher ionization states of different metal species (e.g Steidel et al. 2010;Kacprzak et al. 2015;Chisholm et al. 2018;Ashley et al. 2020), with the hottest phase (> 10 7 ) being accessible in the X-ray emitting gas (e.g.Ptak et al. 1997;Strickland & Heckman 2009;Lopez et al. 2020;Hodges-Kluck et al. 2020).
While there are local starburst systems such as M82 which show evidence of the presence of all of these different phases simultaneously, the situation is more complicated in local dwarf galaxies due to the transient nature of galactic winds and the intrinsically low surface brightness of dwarf galaxy systems (see McQuinn et al. 2019, and references therein for a detailed discussion).Specifically, the presence of a hot wind in nearby dwarfs remains uncertain with only a small fraction of systems showing a hot outflow component (e.g Heckman et al. 1995;Summers et al. 2003Summers et al. , 2004;;Hartwell et al. 2004;Ott et al. 2005;McQuinn et al. 2018), while there is very good observational evidence for warm and cold winds in nearby dwarf galaxy systems (e.g.Meurer 2004;Meurer et al. 2006;Barnes et al. 2011;Kobulnicky & Skillman 2008;Lelli et al. 2014).
Numerical simulations of large cosmological volumes such as Illustris (e.g.Vogelsberger et al. 2014), Magneticum (e.g. Hirschmann et al. 2014), EAGLE (e.g.Schaye et al. 2015), IllustrisTNG (e.g.Pillepich et al. 2018) and SIMBA (e.g.Davé et al. 2019), are not able to directly resolve the processes that drive galactic winds, and therefore rely on phenomenological "sub-grid" prescriptions to implement these winds (see Somerville & Davé 2015;Naab & Ostriker 2017, for reviews).These prescriptions typically involve depositing momentum or energy into the ISM, in a manner that is controlled by adjustable parameters that are tuned to match observations such as the z = 0 stellar mass function (e.g.Baldry et al. 2008Baldry et al. , 2012;;Bernardi et al. 2013;D'Souza et al. 2015), the mass-metallicity relation (e.g.Tremonti et al. 2004), and other observables.While this approach has overall been very successful at reproducing many key observations, it is interesting to note that the required mass loadings are considerably higher than those derived from observations such as the ones referred to earlier (see also Marasco et al. 2022;McQuinn et al. 2019;Xu et al. 2022).However, it is worth noting that the mass loadings derived from observations are highly uncertain, and a true "apples to apples" comparison with simulations (comparing gas outflows with the same phase, measured at the same spatial scale) has not yet been carried out.
It has recently become possible to carry out simulations of sub-galactic portions of the ISM, or even more recently, entire idealized galaxies that can capture the momentum and hot-phase build-up in a resolved Sedov-Taylor phase of individual supernova (SN) remnants (e.g.Kim & Ostriker 2015;Martizzi et al. 2015;Hu et al. 2017;Emerick et al. 2019;Ohlin et al. 2019;Steinwandel et al. 2020;Fujii et al. 2021;Hirai et al. 2021;Grudić et al. 2021).
In these simulations, phenomenological tuning is no longer adopted, and the properties of large scale multiphase winds can be considered emergent predictions following from the physical processes operating on parsec scales.The solar neighborhood environment has been probed extensively in so called "tall box" simulations, which have explored outflow launching mechanisms by supernovae (SNe) either in simulations with fixed SNrate (e.g. de Avillez & Berry 2001;Joung et al. 2009;Hill et al. 2012;Creasey et al. 2013;Fielding et al. 2018) or simulations where the SN-rate follows the star formation rate with a characteristic time-delay (e.g.Girichidis et al. 2016Girichidis et al. , 2018;;Kim & Ostriker 2018;Kannan et al. 2019;Li et al. 2020;Kim et al. 2020a,b;Vijayan et al. 2020;Rathjen et al. 2021;Hu et al. 2021Hu et al. , 2022a;;Kim et al. 2022;Rathjen et al. 2022).
In addition, there have been a number of simulations of wind driving in global galaxies with resolved stellar feedback, but up until recently these were limited to very low mass and low metallicity, isolated dwarf galaxies that resemble the Wolf-Lundmark-Moyette (WLM) dwarf galaxy (e.g.Emerick et al. 2019;Fielding et al. 2017;Hu et al. 2016Hu et al. , 2017;;Hu 2019;Smith et al. 2021;Gutcke et al. 2021;Hislop et al. 2022;Steinwandel et al. 2022), with typical star formation rates between 10 −5 to 10 −3 M yr −1 and star formation rate surface densities between 10 −9 and 10 −6 M kpc −2 yr −1 .While investigations of the outflow driving by SN-feedback have been performed in these kinds of simulations, the fairly low star formation rates under WLM-like conditions and the small system size allow for little variation locally within the galaxy and studies of radial and/or vertical analysis of outflow rates and loading factors are mostly absent from the literature (but see Hu 2019, for an analysis of radially averaged loading factors in WLM-like systems).In addition, the low potential depth of such systems makes it hard to generalize to larger galaxies.
Finally, it is worthwhile to note the immense progress made over the past decade in the regime of "semiresolved" simulations that adopt schemes for resolved SN-injection if the cooling radius is resolved and otherwise input the terminal momentum produced by SNremnants with a given boost factor, often following the previous results of isolated and turbulent blast wave evolution Cioffi et al. (1988); Geen et al. (2016); Kim & Ostriker (2015); Martizzi et al. (2015); Thornton et al. (1998); Walch et al. (2015).Notable examples of these type of simulations are the FIRE-1 (e.g.Hopkins et al. 2014), FIRE-2 (e.g.Hopkins et al. 2018) andFIRE-3 (e.g. Hopkins et al. 2022), the simulations of Kimm et al. (2015), the SMUGGLE model (e.g.Marinacci et al. 2019) as well as the very recent SATIN model (Bieri et al. 2022).In this context it is interesting to point out that these simulations tend to reproduce the rather high mass loading factors in large-scale cosmological simulations.
In this paper we introduce the first results from a simulation of a galaxy that resembles more massive local dwarfs, such as the large Magellanic cloud (LMC) or M33, with resolved stellar feedback, and discuss the outflow launching and driving process in phase space and as a function of radius and vertical height.The structure of this paper is as follows.In Section 2 we will briefly discuss the numerical methods used and describe the initial conditions adopted.In Section 3 we will present the main findings of our study for the time evolution of outflow rates and loading factors, their spatial properties (radial and vertical profiles), as well as the key aspect of this work: the joint PDFs of outflow velocity and sound speed measured at differing heights and radii.In Section 4 we will discuss our results and compare to previous work where appropriate.In Section 5 we will summarize our findings and discuss pathways for future improvements.

Gravity and Hydrodynamics
The simulation discussed in this paper has been carried out with the code p-gadget3 (Springel 2005) with major improvements to the routines for hydrodynamics and stellar feedback included (see Hu et al. 2014Hu et al. , 2016Hu et al. , 2017;;Hu 2019;Steinwandel et al. 2020, for details).The code is a multi-method framework for solving the equations of hydrodynamics and has the option to solve for the fluid flow either with the "pressureenergy" Smoothed Particle Hydrodynamics (SPH Hopkins 2013;Hu et al. 2014) method or the second order meshless finite mass (MFM) method (e.g Lanson & Vila 2008;Gaburov & Nitadori 2011;Hopkins 2015).We use the MFM method as implemented in Steinwandel et al. (2020) which is itself based on the multi-method code gizmo as presented in Hopkins (2015).MFM is realized as a second order Godunov scheme and we solve for fluid fluxes with a Harten-Lax-van-Leer (HLL) Riemann solver with an appropriate reconstruction of the rarefaction wave (HLLC) between tracer particles.This is combined with a second order "kick-drift-kick (KDK)" leap-frog scheme for the time integration.The gravitational forces are computed via the tree particle mesh method (e.g.Barnes & Hut 1986;Springel 2005;Hopkins 2015) and the appropriate corrections for the coupling between MFM and Gadget's gravity solver are applied as outlined in the appendix of Hopkins (2015).

Star formation and stellar feedback
Star formation is implemented in the simulation by sampling the stellar Initial Mass Function (IMF) following the methodology first described in Hu et al. (2016) and Hu et al. (2017) based on a Schmidt-type star formation scaling with a fixed star formation efficiency of 2 per cent per free fall time.In addition, we compute the Jeans mass for all gas cells and convert any gas cell with less than 0.5 M J in the MFM-kernel radius instantly to a star particle in the next time step.The Jeans mass is given via: Additionally, we note that the kernel mass is computed on the average number of neighbors and not the true number of neighbors as computed by the code, to avoid an additional SPH-neighbor search.All star particles above the resolution limit of 4 M are treated as single stars.We follow several feedback channels for the stars with a mass larger than 4 M formed within the simulation, including photoelectric heating and photoionising radiation, as well as the metal enrichment from AGB stellar winds (as a single kernel weighted metal dump at the end of their lifetimes).The core of the stellar feedback scheme is SN-feedback as implemented in Hu et al. (2017) for the SPH version of the code and Steinwandel et al. (2020) for the MFM version of the code, which at the target resolution of our simulation of 4 M allows for a self consistent build-up of the hot phase and the momentum of individual SN-events in a resolved Sedov-Taylor phase (see Hu et al. 2021Hu et al. , 2022a) ) with the deposition of the canonical SN-energy of 10 51 erg into the ambient ISM.Moreover we adopt lifetimes of massive stars based on Georgy et al. (2013), allowing us to self-consistently trace the locations in the ISM where SNe explode.Thus, all the outflow properties in these simulations can be interpreted as being self-consistently launched from the ISM.
Our simulation utilizes a treatment of the far-UV radiation from single stars (above 4 M ) following the implementation of Hu et al. (2017) where we update the local radiation field G 0 during the gravity-tree walk (optically thin approximation).We follow the lifetime functions and effective surface temperatures of Georgy et al. (2013).UV luminosities for each star is taken from the BaSel (Lejeune et al. 1997(Lejeune et al. , 1998;;Westera et al. 2002) stellar library, in the (single bin) energy band 6−13.6 eV, which dominates the photo-electric heating rate in the ISM.We then follow photoelectric heating rates ofBakes & Tielens (1994); Wolfire et al. (2003) and Bergin et al. (2004) with the rate: by multiplying the local radiation field G 0 , the effective attenuation radiation field G eff = G 0 exp(−1.33• 10 −21 cm 2 DN H,tot ), where n denotes the hydrogen number density, D is the dust-to-gas mass ratio, N H,tot is the total (local) hydrogen column density and is the photoelectric heating rate given via: = 0.049 1 + 0.004ψ 0.73 + 0.037(T /10000) 0.7 1 + 2 with ψ = G eff √ T /n e , with the electron number density n e and the gas temperature T .We include a treatment for photo-ionizing radiation since it can have a significant impact on ISM properties and star formation rates (e.g.Hu et al. 2017;Smith et al. 2021;Hislop et al. 2022).Ideally, we would utilize a proper radiative transfer scheme for this purpose, but this is prohibitively computationally expensive.Thus we model the photoionizing feedback via a Strömgren approximation locally around each photo-ionizing source in the simulation, similar to the approach used in Hopkins et al. (2012) and later improved by Hu et al. (2017).
Every star particle with mass m IMF > 8 M is treated as a photo-ionizing source in the presented simulation.All the neighboring gas particles within the Strömgren radius R S of the star: are labeled as photo-ionized, with β as the (electron) recombination coefficient, n the number density and S the photo-ionizing flux.For the gas cells in the Strömgren radius R S we adjust the chemical abundance, destroy all the molecular and neutral hydrogen and transfer it to H + .Furthermore, we heat the surrounding gas within R S to 10 4 K (if it is already warmer than 10 4 K we do not update its temperature), which is the temperature in H II-regions, by providing the necessary energy input from the ionizing source.This procedure is rather simple if R S is a known quantity.In practice, the density structure if the ISM is much more complicated than a constant density background for which R S is well defined and we compute R S in an iterative fashion, by iterating until our ionising photon budget is depleted.
Given this procedure, it is possible that two gas particles are flagged as photo-ionized by two star particles in which case we only adopt the flux of the closer star particle.The current treatment of this procedure is isotropic which leads to a strong mass bias in a scheme.Since we are aware of this we set the maximum of R s to a distance of 50 pc from the host star.This scheme can be improved in the future by healpix tessellating the the ionisation sphere as implemented in Smith et al. (2021).
Our simulation includes a resolved treatment of core collapse supernovae (CCSN) for which the key is ultimately the high mass and spatial resolution with which we run the simulation as pointed out by Steinwandel et al. (2020)for Lagrangian codes one needs a spatial resolution of between 1 and 10 M to resolve the momentum and hot-phase build-up during the Sedov-Taylor phase.Hence our scheme is rather simple and we distribute 10 51 erg with every SNe into the nearest 32 neighbors.We adopt metal enrichment following the yields provided byChieffi & Limongi (2004), who provide data for stellar yields of stars from zero to solar metallicity for progenitors mass from 13 M to 35 M .We interpolate our metal enrichment routines based on this data.When a star explodes in a CCSN the mass is added to the surrounding 32 (one kernel) gas particles with the respective metal mixture that is expected from the above yields.The choice of adopted metal yields is important because it alters the cooling properties of the ambient gas and influences the density and phase structure of the ISM.Finally, we briefly describe how we handle the splitting of particles in our code.At the target resolution we adopt for our simulations one enrichment event can add more mass to the resolution elements that they are originally hold and thus we need to split the new (heavy particles) down to the resolution level of 4 M .Splitting particles is always tricky in Lagrangian codes because it requires some careful handling of the memory available per core to make the procedure memory efficient.We split the gas particles after the enrichment step if the mass exceeds the initial particle mass resolution of the simulation by a factor of two, to avoid numerical artifacts.The split particles inherit all the physical properties that were present in the parent particle including specific internal energy, kinetic energy, metallicity and chemical abundances.The new particle mass is set to one half of the parent particle mass and we assign a random position within one fifth of the parents smoothing length for the new particles position.This process is isotropic within the smoothing kernel of the parent particle.Spawning the new particles can become rather expensive if this has to be carried out repeatedly in the code because if requires a new domain decomposition alongside with a new build of the gravity tree.
Cooling and heating processes are modeled in a similar fashion as in the "SILCC" project (e.g.Walch et al. 2015;Girichidis et al. 2016Girichidis et al. , 2018;;Gatto et al. 2017;Peters et al. 2017) based on a non-equilibrium cooling and heating prescription.The non-equilibrium model is based on the work of Nelson & Langer (1997), Glover & Mac Low (2007a), Glover & Mac Low (2007b) and Glover & Clark (2012).However, we note that some of the rate equations have been updated following Gong et al. (2017).The model tracks six individual species given via molecular hydrogen (H 2 ), ionized hydrogen (H + ), neutral hydrogen (H) as well as carbon-monoxide (CO), ionized carbon (C + ) and oxygen.Additionally, we trace free electrons and assume that all silicon in the simulation is present as Si + .There are several chemical reactions in which these species are involved and we refer the interested reader to table 1 of Micic et al. (2012).The network models molecular hydrogen formation on dust grains for which we solve directly for the non-equilibrium rate equations to obtain x H2 and x H + and obtain x H based on the conservation law: (5) In the network we track n e which we compute via: We adopt non-equilibrium cooling rates that take into account the local density, temperature and chemical abundance.The most dominant cooling processes are fine structure line cooling of atoms and ions (C + , Si + and O), vibration and rotation line cooling due to molecules (H 2 and CO), Lyman-α cooling of hydrogen, collisional dissociation of H 2 , collisional ionisation of hydrogen and recombination of H + (in gas phase and on dust grains).The main heating processes are photoelectric heating from dust grains and polycyclic aromatic hydrocarbonates (PAHs), a constant cosmic ray ionisation rate of 10 −19 s −1 , photo-dissociation of H 2 , UVpumping of H 2 and formation of H 2 .For high temperatures (T > 3 × 10 4 K) we adopt the equilibrium cooling formalism first described by Wiersma et al. (2009) as implemented by Aumer & White (2013) including the elements H, He, C, N, O, Ne, Mg, Si, S, Ca, Fe, and Zn.

Initial Conditions
The initial conditions for this simulation were produced with the method described in Springel et al. (2005), which is typically referred to as "MakeDiskGal" or "MakeGalaxy".As input parameters we choose typical values that are in agreement with observations of the Large Magellanic Cloud (LMC) (e.g.Erkal et al. 2019, based on kinematic arguments) and values previously adopted in numerical simulations to select an LMC-type galaxy (e.g.Hopkins et al. 2018).However, we note that the system is set up to be at the lower limit of halo mass and gas fraction to reduce the computational expense.Hence, we adopt a total dark matter halo mass of 2×10 11 M , a stellar background potential of 2×10 9 M as well as a gas disk with a total mass of 0.5×10 9 M for which we adopt a uniform mass resolution of 4 M , resulting in ∼ 10 8 gas particles (512 3 ) and a mean particle spacing of ∼ 1.0 pc (mean effective spatial resolution) for the gaseous body of the galaxy.For the dark matter halo we adopt a concentration parameter of c = 9 modeled as a Hernquist-profile.We adopt a disc-scale length of 1 kpc and a disc-scale height of 250 pc with an initial metallicity of 0.1 Z .The gravitational softening parameter of the simulations is set to 0.5 pc for gas and stars throughout the total run time of the simulation of 500 Myrs.

ISM and Outflow morphology
In Fig. 1 we show projections of the stellar surface density (left), the gas surface density (center) and a slice of thermal pressure within a thickness of 100 pc (right) for our LMC run in the face-on (top row) and edge-on (center row) view at an evolutionary time of t=400 Myr in a box with a side-length of 10 kpc which is indicated by a white arrow in the bottom of the lower center box.One can clearly see the turbulent structure of the ISM in the density and pressure maps as well as the complex structure of the galaxy in its stellar body.Previous studies of global galaxies with resolved ISM were thus far restricted to WLM-like systems with roughly 10-20 times lower total mass than our LMC-like system and it is worth pointing out some general differences with the much smaller WLM-like systems investigated previously (e.g.Hu et al. 2016Hu et al. , 2017;;Fielding et al. 2017;Hu 2019;Smith et al. 2021;Gutcke et al. 2021;Hislop et al. 2022;Steinwandel et al. 2022).For instance, we note that for the first time we can observe the build-up of more complex structure in the stellar body in a resolved global ISM simulation.This is generally not seen in the lower mass WLM-like systems (or very weakly present depending on the model realization).The face-on view of the gas density reveals large holes in the galaxy, a result of the strong SN-feedback.Additionally, the edgeon view of the gas density reveals large scale outflows with volume filling gas below 0.1 cm −3 and large scale entrained filamentary structure, that is less prominent in the lower mass WLM-runs.The pressure structure of the turbulent ISM is clearly tracing the cavities of the hot bubble-like structures that we find in our simulation and extends out to large radii.
In Fig. 2 we show the radial profiles of the gas and stellar surface density (left hand side) of the young stars with ages less than 40 Myr for t = 300 Myr (solid lines), t = 300 Myr (dashed lines), t = 300 Myr (dash-dotted lines) as well as the star formation rate surface density (right hand side).We note that the system has a central gas surface density of around 3-4 M pc −2 with a stellar surface density that is two orders of magnitude lower.It is interesting to note that star formation is largely restricted to the central part within 3 kpc and is almost flat out to a radius of 2 kpc.This indicates that the wind which we will discuss in the remainder of the paper has a strong impact on the central region of the galaxy from where the wind driving process by SNe originates.

Phase structure of the ISM
In Fig. 3 we show the time-averaged, mass-weighted phase diagrams in density-pressure phase space (left) and in density-temperature phase space (right).The phase diagrams are both evaluated within a slab of 100 pc around the midplane and within a radius of 5 kpc and generally reveal the expected multiphase structure of our simulated ISM.
The warm and cold ISM dominates the total mass budget and covers a wide range of pressure, representing radial variation of the mean heating rate (lower pressure with lower photoelectric heating at outer regions and vice versa).The hot phase is established by SNe.Additionally, we can identify the HII-regions as a horizontal line at T = 10 4 K.We note that there is a significantly higher amount of mass in the diffuse hot phase for our LMC run compared to the smaller WLM-like system of previous studies (e.g.Hu et al. 2016Hu et al. , 2017;;Hu 2019;Steinwandel et al. 2020;Smith et al. 2021;Hislop et al. 2022;Steinwandel et al. 2022).Thus it is particularly interesting to qualitatively compare to the obtained ISM structure of the results of WLM-like systems that also include photo-ionising radiation, but which consistently show less hot and more warm diffuse gas in phase space (e.g.Hu et al. 2017;Smith et al. 2021;Steinwandel et al. 2022).

Definitions for outflow rates and loading factors
We will measure the flow rates at different heights |z|: 0.5 kpc, 1 kpc, 3 kpc, 5 kpc and 10 kpc.Additionally, in some plots we will adopt radial cuts 0.0 < R < 5.0 kpc, 0.0 < R < 1.5 kpc, 1.5 < R < 3.0 kpc and 3.0 < R < 5.0 kpc, where R denotes galactocentric distance R = x 2 + y z .Furthermore, to illustrate the complex structure of the wind it is useful to also measure flow rates in spherical shells, specifically at higher altitude above the midplane.Whenever wind profiles are measured in spherical shells they will be indicated with r, which is defined as r = x 2 + y z + z 2 .Vertical profiles (perpendicular to the midplane) will be indicated as the third component of the gas position vector (z-axis).
We will adopt the following definitions for radial inflow and outflow rates, with the net total flow rates defined via Ṁ = Ṁout − Ṁin and Ė = Ėout − Ėin .Outflow is defined by a positive value of the radial velocity given v • r = v r > 0. We can write down the discrete outflow rates for mass and energy: with γ = 5/3 and dr as the shell thickness in which we measure the outflow rates.We note that for vertical The density-pressure phase space reveals the thermally unstable medium with ∂P/∂n < 0 at a pressure of around 3000 K cm −3 .In the density-temperature phase space we can see that there is clear three phase medium with a hot (SN-heated) phase, a cold (star forming phase) and a warm diffuse phase, that can be split into the warm neutral medium (WNM) and the warm ionized medium (WIM).The horizontal line at 10 4 K represents the photo-ionized gas.
measurement of outflow rates (in plane parallel cylinders at the given heights |z above) we will adopt dr = 0.1 × |z|.For the spherical measurement at r vir , we adopt a spherical shell at 0.1 r vir with a thickness of 0.01 r vir .
It can be useful to investigate these quantities when normalized to a set of reference quantities, defined as the so-called loading factors for mass (η m ), and energy (η e ) that we define via: 1. outflow mass loading factor: η out m = Ṁout /SFR, 2. energy outflow loading factor: In our work we adopt E SN = 10 51 erg, the supernova rate R SN = SFR/(100M ) which are adopted for best comparison and stellar population synthesis modeling with earlier dwarf galaxy studies who adopt these values for normalization (e.g.Hu 2019).We compute the mean star formation rate SFR over a time span of 300 Myr between 200 and 500 Myr of evolution when we investigate the time evolution of loading factors.However, we note that the definition of loading factors in galactocentric, radial or vertical profiles is more difficult to define due to the complicated structure of the wind.
Hence, when we investigate the spatial properties of the wind, we will do this mostly in "outflow rate space" but note that one can always compute a loading factor using the mean star formation rate of SFR ∼ 0.61 M yr −1 and the SN-energy that we provided for the computation of R SN .However, the problem is that specifically at larger heights this approach is not very physical and the correct way to obtain loading factors would be to trace wind streamlines 1 , something we will postpone to future work.Thus the only time that we adopt explicit loading factors in this paper is in Sec.3.7 where we compute the loading factors in bins close to the disc at 0.5 and 1 kpc respectively, where there is still a spatial correlation possible with the star formation activity that is measured "down the barrel" Lastly, we define the Bernoulli velocity as: where c s is the isothermal sound speed, v out is the velocity of the wind, and γ = 5/3 everywhere in the simulation domain.The Bernoulli-velocity represents the terminal velocity that an adiabatic wind would reach in the limit of v B v esc

Time evolution of the outflows
In Fig. 4 we show the time evolution of the star formation rate (top), the mass loading factor η M (center) and the energy loading factor η E (bottom).We compute the star formation rate averaged over four different time intervals by summing up the masses of all stars with ages less than 20, 10, 5 and 1 Myr respectively.While this is merely a sanity check it indicates that the "burstiness" of the star formation history depends on the time interval over which it is averaged.Initially, the star formation rate peaks at a value (first peak) of ∼ 0.1 M /yr,  (black), 10 Myr (blue), 5 Myr (red) and 1 Myr (orange).The center panel shows the time evolution of the mass-loading factor measured at heights of 1 kpc (black), 3 kpc (red), 5 kpc (blue) and 10 kpc (orange).In the bottom panel we show the energy loading as a function of time, again measured at the same four heights above and below the mid plane.The mass loading is normalized to the mean star formation rate of ∼ 0.06 M yr −1 and the energy loading is normalized to the star formation rate per 100 M times 10 51 erg, as the canonical SN-energy, and is thus a measure of how much energy escapes per SN.
which is a remnant from the initial transient phase in this type of isolated galaxy simulation that leads to a small starburst.However, after this transient phase the star formation rate drops by a factor of around 2 and shows a slight increase around 150 Myr followed by a drop towards 250 Myr before a steeper increase around 310 Myrs, which is actually a result of the spatial overlap of several "star clusters" in the lower left quadrant of the galaxy, causing a second (smaller) starburst event.We note that the system has a period of around 150 Myr between peaks in star formation.A similar cycle can be seen in the WLM-like simulations (e.g.Hu et al. 2017;Hislop et al. 2022;Steinwandel et al. 2022) as well as in the Tigress "tall box" simulations (e.g.Kim & Ostriker 2018;Kim et al. 2020b), carried out with the grid code athena (Stone et al. 2008), for their low gas surface density model.
The mass loading and energy loading factors are computed at four different heights: 1 kpc, 3 kpc, 5 kpc, and 10 kpc in slabs above and below the midplane based on the definitions for outflow rates and mass loadings summarized in Section 3.3.Additionally, in Fig. 4 we show results that correspond to a global measurement of the loading factors at 0.1 R vir ∼ 10.0 kpc, similar to the scale often used in studies of gas flows in cosmological simulations (e.g.Pandya et al. 2021).The thickness of the slices at z = 1 kpc, z = 3 kpc, z = 5 kpc and z = 10 kpc correspond to 10 per cent of the height dr = 0.1z, which results in slightly thicker slabs at higher altitude to ensure that the number of particle in the respective slabs is high enough to obtain a meaningful average.Broadly, the peaks in star formation are followed by an increase in the loading factors with a delay time of roughly ∼ 25 Myr at 1kpc, ∼ 30 Myr at 3 kpc, ∼ 35 Myr at 5 kpc and ∼ 50 Myr at 10 kpc.These numbers are consistent between the two peaks in star formation at 50 Myr and 300 Myr that we can observe over the 500 Myrs of evolution of the system.The mass loading is around 1 − 5 at 1 kpc and a factor of a few lower at higher altitude, while the energy loading is generally low with values around a few times of 0.01.

Phase structure of the outflows
A key aspect of this outflow study is to investigate the multiphase nature of galactic outflows.In Fig. 5 we show the joint PDFs of outflow velocity and sound speed weighted by mass outflow rate (left) and energy outflow rate (right) for a radial region between 0.0 < R/kpc < 5.0 at a height of 1 kpc above and below the midplane.Generally, we find that the bulk of the mass is transported by the warm phase of the wind with velocities between 30 km s −1 and 50 km s −1 while for  2021), we find a significant contribution to the energy outflow rate from the warm phase of the wind at least at lower altitude.
In Fig. 6 we show a compilation of the mass outflow rate weighted joint PDFs for different regions above and below the midplane.Different columns indicate different radial regions within the slabs, while different rows indicate different heights above and below the midplane.The PDFs in the first column are measured between 0.0 < R < 1.5 kpc at |z| =0.5 kpc, 1 kpc, 3 kpc, 5 kpc and 10 kpc above and below the midplane in slabs with thickness dr = 0.1z.The other columns show the same respective heights, but are restricted to 1.5 < R < 3.0 kpc (second column), 3.0 < R < 5.0 kpc (third column).Breaking down the distribution shown in Fig. 5 into different regions allows for a more detailed understanding of the complex structure of the outflows in global galaxy models.We highlight some important features here.At a height of 0.5 and 1 kpc we find that the outflow structure has only a very weak radial trend.However, the highest velocities and temperatures in the outflows are reached between 1.5 < R < 3.0 kpc, which is likely related to the fact that the outflowing gas from the center is leaving the galaxy with some finite opening angle.Moreover, the wind carries more mass in the central regions than in the outskirts at lower altitude.However, the effect is modest at 0.5 and 1 kpc and is subsequently inverted at higher altitude, again indicating that the gas accelerated from the center is not ejected perfectly perpendicular to the midplane, most likely due to a combination of the asymmetry of large scale SNbubbles and the large scale shear due to rotation of the galaxy.Considering the distribution in the third column one can observe a quite drastic increase in the high velocity hot component, which can be explained as the hot gas launched from the center that is then adiabatically cooled due to the velocity streamline opening at higher altitude.We note that the area of the rings that we average over is increasing from the center to the outskirts, which might be able to partially explain the inverted radial trend at high altitudes.Using mass fluxes (outflow rate normalized by area) instead of mass outflow rates, we checked that the trend persists.
Furthermore, we observe that there is a gradual loss of the warm component of the wind as we reach higher altitudes and we find that at a height of 10 kpc, only the (cooled) hot component of the wind is left, with most of the mass located in the outskirts of the wind.We find that the low specific energy part of the wind has dropped out.The slow material falls back to the disc.However, it is interesting to note that there is a part of the warm wind that is moving almost as fast as parts of the hot wind, that could potentially be entrained in the hot wind.This would be consistent with the lower sound speed of the wind material at higher altitude, indicating that there is either cooling due to adiabatic expansion of the wind material or because a significant amount of the warm fast part is mixing with the hot wind (see Fielding & Bryan 2022, and references therein for the mixing of cold gas in a hot wind).The exact mixing adiabatic expansion of the wind is an interesting question in itself but would need a more detailed study of the distribution of specific entropy as a function of height.
In Fig. 7 we show the same cuts in height and radius as in Fig. 6 but this time for the energy outflow rate.lar to the case for mass, we find two distinct regions that transport the energy, a low temperature, high velocity part of the wind at a sound-speed of c s ∼ 10 − 15km s −1 with a temperature of around ∼ 10 4 K and a high temperature, high velocity part of the wind that reaches temperatures of up to ∼ 10 7 K and outflow velocities up to a few 100 km s −1 .Qualitatively, we observe similar trends as for the mass outflow rate.At low altitude, just above the galactic midplane, we detect the highest energy flux, while the energy flux drops towards the galactic outskirts at low altitude.In contrast, we observe that the outskirts transport a significantly larger amount of energy than the central part of the galaxy at high altitudes (specifically at 10 kpc height).Furthermore, we again point out that the outflow velocity of the hot increased from 1 kpc to 3 kpc in the outermost ring ranging from 3.0 < R < 5.0 kpc.This is a clear indication that the streamlines of the wind show an open topology as the outflow reaches higher altitudes above and below the midplane (see Sec. 3.6).We note that it is not necessarily surprising that the streamlines of the wind are opening up, but it will have important consequences for the computation of loading factors, which is the main reason why we will focus in the remainder of the paper either on outflow rates directly or loading factors close to the midplane, and postpone a more detailed analysis of loading factors following streamline topology to future work.

Radial and vertical structure of the outflows
In the left panels of Fig. 8 and Fig. 9 we show the radial trends of the mass outflow rate and the energy outflow rate respectively.Again, we plot these at four different heights (1 (black), 3 (red), 5 (blue), and 10 (orange) kpc).We find that at low altitude the mass outflow rate is higher than at high altitude as a function of cylindrical distance R and peaks around a radius of 2 kpc.The energy outflow rate at low altitude is highest in the innermost part (R < 2kpc) of the outflow, after which it sharply declines.However, at higher altitude the outskirts carry more of the total energy in the outflow.This trend indicates that the wind (especially, the hot component) delivers significant energy from inner, lower regions to outer, upper regions.
In the left panels of Fig. 10 and Fig. 11 we show the vertical dependence of the total respective outflow rates for mass and energy for the global disc (0 < R < 5 kpc, olive) as well as three different radial annuli between 0.0-1.5 kpc (cyan), 1.5-3.0(purple) kpc and 3.0-5.0(green) kpc.All of these represent the averaged vertical profiles between 200 and 400 Myr of evolution.We do this to illustrate the different outflow regimes in the center and in the outskirts as a function of height above the midplane, which we could already clearly identify in phase space in Fig. 6 and Fig. 7.We find profiles that decrease with increasing radius at all heights.We note that the sum of the cyan, purple and green lines is the olive line.The variation of the outflow rates and the loading factors gets flatter (but not actually flat) above a height of 5 kpc.The instantaneous profiles are similar to the averaged profiles, modulo small amplitude perturbations that arise from the fluctuations in the star formation rate.This set of figures provides more evidence that the outflows are initially driven from the center and lose mass and energy as a function of increasing height.Most of the mass and the energy are transported to the outskirts of the outflow at higher altitude, while at lower altitude above the midplane the mass budget is comparable in the center and the outskirts.However, the bulk of the energy is initially transported outwards through the center until a height of around 4 kpc, after which the outskirts transport most of the energy.These findings are in good agreement with the detailed phase-space analysis of the outflow in Fig. 6 and Fig. 7.
In the right panels of Fig. 8 and Fig. 9 we show the radial trends of the mass outflow rate and the energy outflow rate respectively separated by the warm (T < 5 × 10 5 , dashed lines) and hot phase respectively (T > 5 × 10 5 , dotted line) for four different height above the midplane at 1 (black), 3 (red), 5 (blue) and 10 kpc (orange).The mass outflow rate as function of the radius R is dominated at all heights that we investigated by the warm wind and most of the mass (∼ 90 per cent) is transported by this component while the hot component transport ∼ 10 per cent of the total mass.However, at low altitude most of the mass in the hot wind is transported through the central "chimney", while at higher altitude most of the hot gas in the wind is located in the outskirts of the galaxy, indicating that the hot wind is leaving trough an angle and is mostly generated in the central starburst region, while the warm wind shows almost constant outflow rates as function of radius.
In the right panels of Fig. 10 and Fig. 11 we show the vertical trends of the mass outflow rate and the energy outflow rate respectively separated by the warm (T < 5 × 10 5 , dashed lines) and hot phase respectively (T > 5×10 5 , dotted line) for four different radial annuli along the midplane given via 0.0 < R < 1.5 (cyan), 1.5 < R <  (200 Myrs; radial profiles for the total mass outflow rate at different heights of 1 kpc (black), 3 kpc (red), 5 kpc (blue) and 10 kpc (orange) on the left and the mass outflow rate split in warm (T < 5 × 10 5 K, dashed lines) and hot (T > 5 × 10 5 K, dotted lines) on the right.Generally, the mass outflow rate is higher at lower altitude and exhibits a peak at around R =1.5 to 2 kpc at low altitude.However at 10 kpc we find a gradually increasing trend outwards.If we look a the split in the warm and the hot component on the right it becomes clear that at all radii and at the four different heights at which we computed the radial profiles, we find that the mass transport is dominated by the warm component, where we find that at low altitude (1kpc) roughly 10 per cent of the total outflow rate is in the hot phase, while it appears that there is slightly more mass in the hot at the higher altitudes up to a maximum of around 20 per cent.(200 Myrs; radial profiles for the total energy outflow rate at different heights of 1 kpc (black), 3 kpc (red), 5 kpc (blue) and 10 kpc (orange) on the left and the energy outflow rate split in warm (T < 5 × 10 5 K, dashed lines) and hot (T > 5 × 10 5 K, dotted lines) on the right.The energy outflow rate behaves similarly apart from the fact that there is a sharp drop off beyond a radius of 2 kpc at a height of 1 kpc.In combination with the radial trends of the mass outflow rate presented in Fig 8, these results indicate that the wind in the low altitude, inner regions has moved outward as it propagates upward (especially, energy-carrying hot component); i.e., the wind has an open streamline structure.The split in warm and hot wind on the right reveals that the energy transport is dominated by the hot phase for 1 kpc, 3 kpc and 5 kpc.At an altitude of 10 kpc the warm component is dominating the energetics of the wind. Hwever, this is likely due to the fact that the formerly hot component of the wind is cooled adiabatically expanding as it is expanding. .Vertical profiles for the total mass outflow rate (left), measured in different radial regions given via 0.0 < R < 1.5 kpc (cyan), 1.5 < R < 3.0 kpc (purple), 3.0 < R < 5.0 kpc (green) and 0.0 < R < 5.0 kpc (olive).All of these profiles are averaged over 200 Myrs from 200 to 400 Myrs. A low altitude out to around 2 kpc height the wind is transporting a similar amount of mass in the inner and outer parts of the galaxy.Beyond 2 kpc we find that most of the mass is transported in the outer parts of the wind, indicating an opening velocity streamline structure.If this is split into warm (T < 5 × 10 5 K, dashed lines) and hot (T > 5 × 10 5 K, dotted lines, right hand side), it becomes clear the warm wind is completely dominating the mass transport out to a height of 10 kpc. .Vertical profiles for the total energy outflow rate (left), measured in different radial regions given via 0.0 < R < 1.5 kpc (cyan), 1.5 < R < 3.0 kpc (purple), 3.0 < R < 5.0 kpc (green) and 0.0 < R < 5.0 kpc (olive).All of these profiles are averaged over 200 Myrs, between 200 and 400 Myrs of evolution.At low altitude (1 kpc) most of the energy is transported in the central part of the wind, while at higher altitude we find that there is more energy in the outer parts of the wind.Separating this by warm (T < 5 × 10 5 K, dashed lines) and hot (T > 5 × 10 5 K (right) gas indicates that at low altitude most of the energy is transported by the hot wind out to 2 kpc, in the innermost 2 rings we defined, out to a height of around 7 kpc when the warm wind starts to dominate the energy transport.In the outermost ring the warm wind dominates out to 2 kpc when the hot wind starts to dominate out to 7 kpc after which the warm dominates again.This turnaround in the outermost ring is likely the contribution from the high angle (hot) material from the central wind.Mean values for the radial velocity (top left, vr > 0), vertical outflow velocity (top right), the Bernoulli velocity (bottom left), and the opening angle of the wind (defined as 360/2π × arccos(vz/vr)) within different radial regions between 0.0 < R < 1.5 kpc (black), 1.5 < R < 3.0 kpc (red), 3.0 < R < 5.0 kpc (blue) and 0.0 < R < 5.0 kpc (orange).We note that these results are mass weighted, while our PDFs in Fig. 6 and Fig. 7 are mass outflow rate or energy outflow rate weighted respectively.To demonstrate the effect of weighting either by mass or energy outflow rate we add the mass and energy flux weighted velocities for 0.0 < R < 5.0 kpc as dashed orange (mass outflow rate) and dotted orange (energy outflow rate) lines to the top left panel.
3.0 (purple), 3.0 < R < 5.0 (green) and 0.0 < R < 5.0 kpc (olive).In all radial bins and at all height the mass outflow is dominated by the warm wind, which is especially true at low altitude out to 2 kpc after which the mass fraction between hot and warm remains close to 10 per cent.At low altitude out to 500 pc the sharp increase in mass outflow rate in the warm is likely related to the fact that we pick up a lot of warm material that belongs to the disc and classify it as outflow.The energy outflow rate on the other hand is dominated by the warm wind out to a height of 7 kpc in all radial bins.This is likely due to the fact that the hot wind is adiabatically cooled with increase in height and the formerly hot wind is classified as warm wind above this specific height.In future studies this has to be investigated more closely by the profiles for specific entropy to disentangle the origin of the reduction of temperature of the hot out to higher heights which is beyond the scope of this paper which has the main focus on the phase separation of outflows.
In Fig. 12 we show the time averaged (over 200 Myr) radial outflow velocity (top left), the time averaged vertical outflow velocity (top right), the time averaged Bernoulli velocity (bottom left) and the time averaged opening angle (bottom right) of the wind as a function of height above the midplane.We find that all averaged velocity profiles increase steeply until a height of around 2.5 to 3.5 kpc up to values of around ∼100 km s −1 , where the Bernoulli velocity reaches slightly higher values than the outflow velocity (see the definition of v B ). Generally, these results are very comparable to 13. Top: Edge-on slice at y = 0 with a thickness of 50 pc of hydrogen number density at t = 330 Myr, demonstrating the large scale outflow from -20 kpc to 20 kpc in x and z direction indicating an opening streamline geometry.Bottom: Zoom on the innermost disc (-5 to 5 kpc width and -1.5 kpc to 1.5 kpc height) to capture parts of the smaller scale (turbulent) fountain flow.The vectors show the directionality (not magnitude) of the velocity field and clearly identify that the velocity field is opening up as a function of height.
those of previous studies such as Kim et al. (2020a).The solid lines in the panels in Fig. 12 for radial, vertical and Bernoulli velocities are mass weighted.We state this specifically to avoid confusion with Fig. 6 and Fig. 7 which are the mass outflow rate and energy outflow rate weighted distributions (they are shifted by v r ).Hence, if we had computed the mean outflow velocity using the weighting from Fig. 6 and Fig. 7 respectively we would obtain a velocity profile with a higher mean velocity.We demonstrate that effect on the time averaged radial outflow velocity for 0 < R < 5.0 kpc with the dashed (mass flux weighted) and the dashed-dotted line (energy flux weighted).
Finally, we note that the shape of the opening angle in the bottom right panel of Fig. 12 is increasing as a function of height, peaks around 2.5 to 3.0 kpc and declines outwards to a height of 10.0 kpc.We compute the opening angle α via: Qualitatively, this shape can be understood as the different components of the wind.To visualize this we plot the streamlines of the wind in Fig. 13 at t = 330 Myrs corresponding to the peak in outflow rate following the peak in star formation at t = 300 Myrs.In the top panel we show the large scale gas outflow out to 20 kpc height above and below the mid plane.The white arrows show the direction of the gas itself.We generally can say that, in the central part, the opening angles are small because the wind is accelerated almost vertically.At larger radii, the opening angle gets larger overall because the bulk of the wind is launched from the innermost 3.0 kpc, the region of the largest star formation rate.The streamlines are then bent as the outflow widens and the opening angle increases at larger radii.The turnover towards decreasing opening angles at larger heights occurs because material with large opening angles will simply not reach these heights at the given radial bin and subsequently "rain down" on the outer galactic disc.In the bottom panel of Fig 13 we visualize the small scale turbulent flow closer to the mid plane indicating the fountain flow on smaller scales.

Loading factor dependence
In Fig. 14 we show the dependence of the mass loading factor η M (top) and the energy loading factor η E (center) for a height of 0.5 kpc (left) and a height of 1.0 kpc (right), on the star formation surface rate density.The data in these plots are computed on a grid with a cell size of 2.5 × 2.5 kpc 2 resulting in 16 cells across the innermost 5 kpc of the galaxy.The different colors represent the warm (blue) and hot (red) loading factors, where the warm and hot material are separated at 5 × 10 5 K.The data in these plots is collected over 20 snapshots from 200 Myrs to 400 Myrs and the star formation rate surface density is computed as a 40 Myr time average.For the mass outflow rate we find that at low altitude there is a weak correlation between the star formation rate surface density and the outflow rate, while there is almost no correlation with the hot mass loading.We provide the fits for the loading factors obtained at heights of 0.5 and 1.0 kpc, for the mass and energy in the warm and hot wind respectively.The scalings we find are as follows, where we give Σ sfr,40Myr in the following equations in units of M yr log η E,warm (0.5kpc) = −0.14log Σ sfr,40Myr M yr −1 kpc −2 −2.14 log η E,hot (0.5kpc) = +0.13log Σ sfr,40Myr M yr −1 kpc −2 − 1.20 (16) log η E,warm (1kpc) = −0.20 log Σ sfr,40Myr M yr −1 kpc −2 − 2.68 (17) log η E,hot (1kpc) = +0.023log Σ sfr,40Myr M yr −1 kpc −2 − 1.82 (18) Furthermore, we directly compare to the fits found in the tigress simulation suite based on the results of Kim et al. (2020b).We note that we generally find good agreement with the slopes of the scalings from tigress for the warm wind, while the normalization in our simulation is lower by 0.5-1 dex.Additionally, there are some notable differences for the hot wind in our simulations.The reasons for this will be discussed below.

DISCUSSION
In this section we discuss the implications of the results presented in Section 3, compare to previous work, and describe limitations and pathways for future and follow up work.

The multiphase ISM
Fig. 1 shows a turbulent galactic ISM with large SNdriven bubbles all over the midplane.The strong stellar feedback is generating large holes in the ISM with diameters of up to one kpc.The observed structure of the multiphase ISM is similar to that of previous studies of WLM-like dwarfs such as presented in Hu et al. (2016Hu et al. ( , 2017)); Hu (2019), Lahén et al. (2019aLahén et al. ( ,b, 2022)), Smith et al. (2021), Gutcke et al. (2021) and Hislop et al. (2022).From these models the simulations of Hu et al. (2017) labeled as SN-PI-PE have the most comparable physical setup.However, our dwarf is an LMC-like analogue that is 10 times more massive (gas mass and halo mass)than the WLM-like simulation of Hu et al. (2017).Generally, we find that the phase structure of the ISM in density-temperature phase space as presented in the bottom panel of Fig. 3 looks very similar to the results of Hu et al. (2017) with the difference that we find more mass in the hot phase in our LMC-run.Hence, our density-temperature phase diagram shows more similarities to the phase diagram of Lahén et al. (2019a) who simulated a merger of two WLM-like dwarfs that shows a peak of star formation with similar values to those we find in our LMC-like run of around 0.1 M yr −1 based on our Fig. 4. The same can be said when comparing their density-pressure phase space diagram with ours, which shows a similar thermally unstable region around a density of 0.1 cm −3 and pressure of 3000 K cm −3 as in the merger simulations of Lahén et al. (2019a) and Lahén et al. (2019b).The studies of Smith et al. (2021) and Hislop et al. (2022) put the focus on a more compact version of the WLM-like dwarf galaxy with a factor of 20 lower average star formation rate than our system.While Hislop et al. (2022) did not directly look into the ISM-phase structure, the initial conditions and the code version used are the same as in the compact WLM version of Hu et al. (2017) and similar results are expected.The phase-structure of Smith et al. (2021) looks slightly different from the one we find, but they also use a different cooling and chemistry network based on grackle (Smith et al. 2017).The same can be said for the simulations published by Gutcke et al. (2021), who use a tabulated cooling function and a density threshold for star formation.We note that a more comprehensive comparison between different hydrodynamical methods can be found in Hu et al. (2022b) for the WLM regime.

The structure of the outflow
The central result of this study is the morphology and phase structure of the outflow driven by SN-feedback which we report for the first time in a simulation of an LMC-like analogue at a mass resolution of 4 M , which allows us to self-consistently resolve SN-feedback within the multiphase ISM.Some of the first simulations of global galaxy simulations with improved treatment of the ISM-physics were carried out by Hopkins et al. (2012), who simulated four galaxies in the mass range of a Milky Way (MW), an SMC, a gas rich Sbc and a thicker disc that represented higher redshift conditions.The run we present here is in between the runs SMC and Sbc of Hopkins et al. (2012), where we note that the SMC model of Hopkins et al. (2012) is a factor of 1.5 more massive (in gas) than our LMC run based on the gas mass adopted, while our LMC is a scaled up version of our previous WLM runs (Hu et al. 2017;Steinwandel et al. 2022).Additionally, Hopkins et al. (2012) measured mass loading factors based on a positive Bernoulli parameter which they define as b = (v 2 out + 3c s − v esc ) and compared to a definition where the outflow is comprised out of gas with v z > 100 km s −1 and |z| > 500 pc, where they compute a mean by M out / M , and report similar number for these procedures.Our definition is based on fluxes in plane parallel slabs at different height.We find a similar value of around 1 to 5 for eta eta at a height of 500 pc for the time evolution of the mass loading factors as shown in the center panel of Fig. 4 when directly compared to Fig. 6 of Hopkins et al. (2012) for their SMC and Sbc runs.Our mass loading factors decrease by around 1.5 dex at higher altitude.
More recently, Pandya et al. (2021) performed a detailed study of mass, metal, momentum and energy loading in the FIRE-2 cosmological zoom-in simulations.The most comparable halos to our simulation are the FIRE-2 m11 halos which show higher mass loading factors compared to the ones that we report based on our results in the center panel of Fig. 4 within 0.1 R vir .The energy loading reported in Pandya et al. (2021) is higher than our time averaged global values that we find based on the results of the bottom panel of Fig. 4, where we find an energy loading of 0.01 to 0.07 at 1 kpc height, but if we measure at R vir where Pandya et al. (2021) measured their energy loadings we find a value that is consistently around 0.01 while Pandya et al. (2021) report values from 0.05 to 0.3 for the systems that come closest to our simulations (m11a, m11b, m11q, m11c, m11v, m11f).
One of the central parts of our study is the detailed analysis of the 2d joint PDFs of outflow velocity and sound speed for the mass flux (Fig. 6) and the energy flux (Fig. 7).These have been investigated previously in a number of studies including Fielding et al. (2018), Kim et al. (2020a), Schneider et al. (2020) and Steinwandel et al. (2022) with the general conclusion that the bulk of the mass is transported by the warm phase while most of the energy is transported in the hot phase of the wind.However, Pandya et al. (2021) pointed out that the hot outflow loading dominates for M > 10 8 M Fielding et al. ( 2018) and Kim et al. (2020a) studied these PDFs for vertically stratified galactic discs and we note that we generally find good agreement between the shape of their 2d PDFs and ours, when compared at the appropriate heights (200pc and 540 pc for Fielding et al. (2018) and at the scale height H, and 2H as well as at 500 pc and 1000 pc for Kim et al. (2020a)).These heights are in the range of the first two rows of our Fig.6 and Fig. 7 at 500 pc and 1000 pc, for which we find quite good agreement on the mass outflow rate where we can clearly see that the bulk of the mass is transported by the warm wind.However, the energy flux shows a more bimodal distribution, with the warm component gradually vanishing at larger heights.Ultimately, this indicates two concepts.First, the slower part of the warm wind com-ponent is falling back to the disc (the left part of the phase diagram vanishes at larger heights).Second, the faster part of the warm wind component (v out > 50 km s −1 ) is lost at higher heights, most likely due to mixing with the hot wind or its adiabatic expansion which would both be consistent with the lower temperature of the hot wind at larger height.In addition, our simulation allows for a more detailed investigation of the shape of the wind by quantifying the radial structure as we described in Section 3. We note that the comparison with the results of Pandya et al. (2021) based on the FIRE-2 simulations (Hopkins et al. 2018) is more complex since they measured the outflow rates in shells between 0.1 and 0.2 of R vir .Qualitatively, the simulations indicate a similar phase structure.
Finally, it is worth discussing the complex structure of the wind in greater detail, for which we find multiple indicators that there has to be an opening velocity stream-line trajectory.The first indicator for this can be seen based on the radial outflow rate profiles of Fig. 8 and Fig. 9 for different heights above the midplane, which shows that the energy outflow rate is decreasing at low height as a function of radius and slightly increasing at the largest height of 10 kpc.This is evidence that there is a redistribution of energy from the center at low heights to the outskirts at large heights.If we consider vertical profiles of the outflow velocities that are of significance for this study we can try to obtain a first order estimate of the opening angle of the wind based on the misalignment of the radial velocity and its vertical component for which we show radial profiles in Fig. 12.These can be used to obtain the opening angle which we measure at different radii as a function of height above the mid plane, where we find the largest opening angle in the outermost regions above the disc and the lowest opening angles in the central parts, as expected.
The shape of the opening angle is quite complex, but one can straightforwardly show that a similar shape can be obtained by assuming a ballistic wind with uniform streamlines and constant angle dependent injection velocity.The idea is then that, in the inner parts of the galaxy at low altitude one is simply measuring a wind component perpendicular to the disc.Since the wind is opening up as a function of height we measure first a larger opening angle until around a height of 3 kpc, after which we find a decline again, simply because most of the bent field lines will not reach that altitude, at least not within the innermost radial region above the disc.This then explains the increase in the opening angle in the different radial regions as we measure more and more of the material that is driven out of the center.
We also find that the wind shows similar scalings for the dependence of the loading factors as a function of the star formation rate surface density in comparison to previous work.The derived scalings for the warm wind are in good agreement with the results of Kim et al. (2020b).For the hot wind we find similar slopes for the mass and the energy loading but the loading factors are lower than in the simulations of Kim et al. (2020b).That being said, the mass loading factors at low altitude are comparable for the warm wind and only slightly lower for the hot wind than the ones reported by Kim et al. (2020b).The energy loading at low altitude is more dominated for our set of simulations by the warm component when compared to Kim et al. (2020b) in the sense that we have similar energy loadings in warm and hot phase with higher laodings in the warm than Kim et al. (2020b) and lower loadings in the hot.At higher altitude we find an offset by half a dex for the mass loading (warm and hot), while the energy loading is consis-tently lower in the hot than in Kim et al. (2020b) but the warm loadings are comparable.One possible explanation for this could be the suppression of clustering of the stars and subsequently the feedback (e.g., Fielding et al. 2017Fielding et al. , 2018)), because our simulations include photo-ionizing feedback that has been shown in previous simulations to decrease the clustering of the SNe (Hu et al. 2017;Smith et al. 2021) and is a component that was not active in the simulations of Kim et al. (2020b) Additionally, we note that the lower energy loading factors are generally also found in the lower mass systems we simulated such as WLM and LeoP which will be presented with a detailed discussion on the wind phasestructure in a future study.Specifically, when considering the results of Steinwandel et al. (2022), it is quite clear that there is almost no hot wind in the WLMfiducial run in the classical sense and a hot wind was only able to be established by the inclusion of runaway stars.Alternatively, one can increase the gas surface density of the simulations.
We note that Andersson et al. (2022) find a weaker impact of runaways stars on the phase structure of galactic outflows in WLM-like simulations, but they also simulated a higher column density (shorter disc scale length) version of the dwarf that is presented in Steinwandel et al. (2022) with a different physics implementation.
Recently, Rathjen et al. (2022) investigated the phase structure of the winds in the SILCC simulations.When only SNe, stellar winds and photo-ionizing radiation are included they find similar results for the PDFs for mass outflow rate and energy outflow rate in terms of outflow velocities (bulk mass outflow between 10 and 100 km s −1 ) and thermal state of the gas (warm gas with c s ∼ 10 km s −1 ) comparable to what we find.While we do not include a proper treatment for stellar winds, we adopt a treatment for the photo-ionizing radiation, which makes our runs quite comparable to theirs considering the underlying gas cooling is almost identical between SILCC and our simulations, although they adopt higher metallicity.Its unlikely that stellar winds would change this picture strongly given recent work on stellar wind bubbles (e.g.Lancaster et al. 2021a,b,c).Additionally, a Lagrangian code such as ours is not necessarily best suited for a stellar wind feedback implementation due to the lack of resolution elements with which one can resolve the bubble.Resolving the thermal mixing layer of the bubble which is ultimately driving its evolution requires much higher resolution.Hence, any stellar wind treatment at our resolution has to be labeled as "subgrid".Furthermore, Rathjen et al. (2022) investigate the impact of cosmic ray (CR) diffusion driven winds on the phase structure, which leads to an additional cold phase wind launched by CRs above the midplane.Since, we do not adopt this mechanism it is obviously absent from our phase structure.That being said it is not true for our simulations that we do not obtain any cold material in our phase structure below a sound speed of 10 km s −1 , but it is difficult to properly compare since they exclusively select the snapshots which have mass loading higher than one in their analysis, while we average over a full wind driving cycle.Finally, we note that they simulated much higher column densities in their main paper than we did although they show some comparisons to their lower column density runs in the appendix (which are still higher column density than our LMC analogue by at least half dex).

CONCLUSIONS
Finally, we summarize the main findings of this work.We present initial results from an LMC-like simulation with solar mass and sub-parsec resolution with a special focus on the supernova-driven galactic outflows.This is the largest galaxy simulated with a resolved ISM to date. 1.Our system shows a star formation rate that varies relatively little over the evolution time, with a mean star formation rate of around ∼ 0.05 M yr −1 and a cycle of roughly 150 Myrs.The values for gas surface density and star formation rate surface density we put forward in Fig. 2 indicate that the galaxy is forming stars in agreement with the observations of Bigiel et al. (2010) for the dwarf galaxy regime, corresponding to a gas depletion time scale of 10 Gyr.
2. We measured the time evolution of loading factors for mass and energy.At low altitude (1 kpc) we find a mass loading around ∼ 1.0 while at higher altitudes we find a roughly constant global value of around 0.1.The global averages for the energy loading are roughly constant with a slight decrease from 2 × 10 −2 to 9 × 10 −3 as we move outwards from 1 kpc to 10 kpc.This indicated that the wind at higher altitude is dominated by the hot wind.
3. We discussed the mass and energy flux weighted joint PDFs of outflow velocity and sound speed in different radial regions and different heights.
We see a bimodal wind structure, where both the warm and hot wind transport significant amounts of mass and energy.These results are qualitatively similar to the findings of Kim et al. (2020a) as well as the simulations of Fielding et al. (2018) but differ in terms of the total mass and energy loading budget as in our simulations hot and warm mass and energy loading are more similar to one another than in Kim et al. (2020a) and Fielding et al. (2018).
4. The wind shows an opening streamline structure with an average opening angle of around 30 degrees that has a peak at a height of around 3 kpc.The structure of the streamlines is in rough agreement with a (uniform) ballistic wind.However, a more detailed analytical model for the streamline structure needs to be developed to better understand the complex structure of the opening angle.
5. The mass loading of the warm wind is decreasing with the star formation rate surface density Σ SFR with a slope of around -0.5 at 500 pc and 1000 pc height.The hot wind shows a scaling with Σ SFR with a slope of around +0.01 at 500 pc and -0.06 at 1000 pc.The slopes of the warm and hot wind scaling are in good agreement with the results of Kim et al. (2020a) but both wind components show lower loading factors by roughly half a dex.The energy loading shows asymmetric behavior for the hot and the warm wind with slopes of -0.14 for the warm and +0.13 for the hot wind at 500 pc height with similar trends at 1000pc height for the warm wind, but shallower for the hot wind (slope of +0.023).The hot wind has lower energy loading and the warm wind has higher energy loading than in Kim et al. (2020b) at similar star formation rate surface density.There are a number of reasons for why exactly that could be the case, including the different gravitational potential between our runs and tigress or the different geometry adopted.Hence, this is not necessarily to be interpreted as an apples to apples comparison to tigress and could be explained by testing against pressure regulated feedback modulated theory of Ostriker & Kim (2022) which we postpone to a future study.

Figure 1 .Figure 2 .
Figure1.Face-on (top) and edge-on (bottom) projections at t=400 Myr for stellar surface density (left), gas surface density (center) as well as a slice within 100 pc of the thermal pressure (right).The top center panel clearly shows the extended SN-bubbles all across the galaxy.The edge-on view reveals a large scale outflow, easily reaching out to a height of 5 kpc with volume filling gas of density around 0.1 cm −3 .The face-on projection is restricted to 500 pc above and below the mid-plane.The edge-on view is a full projection over 10 kpc along the line of sight in order to visualize the low density gas in the outflows.

Figure 3 .
Figure3.Mass weighted structure of our multiphase ISM in density-pressure (top) and density-temperature (bottom) phase space.The density-pressure phase space reveals the thermally unstable medium with ∂P/∂n < 0 at a pressure of around 3000 K cm −3 .In the density-temperature phase space we can see that there is clear three phase medium with a hot (SN-heated) phase, a cold (star forming phase) and a warm diffuse phase, that can be split into the warm neutral medium (WNM) and the warm ionized medium (WIM).The horizontal line at 10 4 K represents the photo-ionized gas.

Figure 4 .
Figure 4.The top panel shows the time evolution of the star formation rate, averaged over different timescales of 20 Myr(black), 10 Myr (blue), 5 Myr (red) and 1 Myr (orange).The center panel shows the time evolution of the mass-loading factor measured at heights of 1 kpc (black), 3 kpc (red), 5 kpc (blue) and 10 kpc (orange).In the bottom panel we show the energy loading as a function of time, again measured at the same four heights above and below the mid plane.The mass loading is normalized to the mean star formation rate of ∼ 0.06 M yr −1 and the energy loading is normalized to the star formation rate per 100 M times 10 51 erg, as the canonical SN-energy, and is thus a measure of how much energy escapes per SN.

Figure 5 .
Figure5.2d PDFs of outflow velocity and sound speed color coded by mass outflow rate (left) and energy outflow rate (right), at a height of 1 kpc for the radial region between 0.0 < R < 5.0.We note that the left panel is the sum of the three panels of the second row in Fig.6and the right panel is the sum of the three panels of the second row in Fig.7.The grey lines in the PDFs define contours based on equal Bernoulli-velocity as defined by equation eq. 9.All the PDFs presented in this paper are time averaged between 200 and 400 Myrs of evolution.theenergy we find a fast hot component that carries at least 50 percent of the energy of the wind.While this is similar to previous results as shown for example byFielding et al. (2018),Kim et al. (2020a), orPandya et al. (2021), we find a significant contribution to the energy outflow rate from the warm phase of the wind at least at lower altitude.In Fig.6we show a compilation of the mass outflow rate weighted joint PDFs for different regions above and below the midplane.Different columns indicate different radial regions within the slabs, while different rows indicate different heights above and below the midplane.The PDFs in the first column are measured between 0.0 < R < 1.5 kpc at |z| =0.5 kpc, 1 kpc, 3 kpc, 5 kpc and 10 kpc above and below the midplane in slabs with thickness dr = 0.1z.The other columns show the same respective heights, but are restricted to 1.5 < R < 3.0 kpc (second column), 3.0 < R < 5.0 kpc (third column).Breaking down the distribution shown in Fig.5into different regions allows for a more detailed understanding of the complex structure of the outflows in global galaxy models.We highlight some important features here.At a height of 0.5 and 1 kpc we find that the outflow structure has only a very weak radial trend.However, the highest velocities and temperatures in the outflows are reached between 1.5 < R < 3.0 kpc, which is likely related to the fact that the outflowing gas from the center is leaving the galaxy with some finite opening angle.Moreover, the wind carries more mass in the central regions than in the outskirts at lower altitude.However, the effect is modest at 0.5 and 1 kpc and is subsequently inverted at higher altitude, again indicating

Figure 6 .
Figure6.Compilation of the 2d joint PDFs for mass outflow rate (mass loading) for different regions in the galaxy.The columns from left to right show the PDFs in different radial regions 0.0 < R < 1.5, 1.5 < R < 3.0 and 3.0 < R < 5.0.The rows from top to bottom show the joint PDFs at heights of 0.5, 1, 3, 5 and 10 kpc in slabs with a thickness of 50, 100, 300, 500 and 1000 pc, between 0.0 < R < 5.0 kpc.The grey dashed lines represent surfaces of equal Bernoulli-velocity following eq.9.All PDFs are averaged between 200 and 400 Myr of evolution.

Figure 8 .
Figure8.Time averaged(200 Myrs; radial profiles for the total mass outflow rate at different heights of 1 kpc (black), 3 kpc (red), 5 kpc (blue) and 10 kpc (orange) on the left and the mass outflow rate split in warm (T < 5 × 10 5 K, dashed lines) and hot (T > 5 × 10 5 K, dotted lines) on the right.Generally, the mass outflow rate is higher at lower altitude and exhibits a peak at around R =1.5 to 2 kpc at low altitude.However at 10 kpc we find a gradually increasing trend outwards.If we look a the split in the warm and the hot component on the right it becomes clear that at all radii and at the four different heights at which we computed the radial profiles, we find that the mass transport is dominated by the warm component, where we find that at low altitude (1kpc) roughly 10 per cent of the total outflow rate is in the hot phase, while it appears that there is slightly more mass in the hot at the higher altitudes up to a maximum of around 20 per cent.

Figure 9 .
Figure9.Time averaged(200 Myrs; radial profiles for the total energy outflow rate at different heights of 1 kpc (black), 3 kpc (red), 5 kpc (blue) and 10 kpc (orange) on the left and the energy outflow rate split in warm (T < 5 × 10 5 K, dashed lines) and hot (T > 5 × 10 5 K, dotted lines) on the right.The energy outflow rate behaves similarly apart from the fact that there is a sharp drop off beyond a radius of 2 kpc at a height of 1 kpc.In combination with the radial trends of the mass outflow rate presented in Fig8, these results indicate that the wind in the low altitude, inner regions has moved outward as it propagates upward (especially, energy-carrying hot component); i.e., the wind has an open streamline structure.The split in warm and hot wind on the right reveals that the energy transport is dominated by the hot phase for 1 kpc, 3 kpc and 5 kpc.At an altitude of 10 kpc the warm component is dominating the energetics of the wind.However, this is likely due to the fact that the formerly hot component of the wind is cooled adiabatically expanding as it is expanding.
Figure10.Vertical profiles for the total mass outflow rate (left), measured in different radial regions given via 0.0 < R < 1.5 kpc (cyan), 1.5 < R < 3.0 kpc (purple), 3.0 < R < 5.0 kpc (green) and 0.0 < R < 5.0 kpc (olive).All of these profiles are averaged over 200 Myrs from 200 to 400 Myrs.At low altitude out to around 2 kpc height the wind is transporting a similar amount of mass in the inner and outer parts of the galaxy.Beyond 2 kpc we find that most of the mass is transported in the outer parts of the wind, indicating an opening velocity streamline structure.If this is split into warm (T < 5 × 10 5 K, dashed lines) and hot (T > 5 × 10 5 K, dotted lines, right hand side), it becomes clear the warm wind is completely dominating the mass transport out to a height of 10 kpc.
Figure11.Vertical profiles for the total energy outflow rate (left), measured in different radial regions given via 0.0 < R < 1.5 kpc (cyan), 1.5 < R < 3.0 kpc (purple), 3.0 < R < 5.0 kpc (green) and 0.0 < R < 5.0 kpc (olive).All of these profiles are averaged over 200 Myrs, between 200 and 400 Myrs of evolution.At low altitude (1 kpc) most of the energy is transported in the central part of the wind, while at higher altitude we find that there is more energy in the outer parts of the wind.Separating this by warm (T < 5 × 10 5 K, dashed lines) and hot (T > 5 × 10 5 K (right) gas indicates that at low altitude most of the energy is transported by the hot wind out to 2 kpc, in the innermost 2 rings we defined, out to a height of around 7 kpc when the warm wind starts to dominate the energy transport.In the outermost ring the warm wind dominates out to 2 kpc when the hot wind starts to dominate out to 7 kpc after which the warm dominates again.This turnaround in the outermost ring is likely the contribution from the high angle (hot) material from the central wind.
Figure12.Mean values for the radial velocity (top left, vr > 0), vertical outflow velocity (top right), the Bernoulli velocity (bottom left), and the opening angle of the wind (defined as 360/2π × arccos(vz/vr)) within different radial regions between 0.0 < R < 1.5 kpc (black), 1.5 < R < 3.0 kpc (red), 3.0 < R < 5.0 kpc (blue) and 0.0 < R < 5.0 kpc (orange).We note that these results are mass weighted, while our PDFs in Fig.6and Fig.7are mass outflow rate or energy outflow rate weighted respectively.To demonstrate the effect of weighting either by mass or energy outflow rate we add the mass and energy flux weighted velocities for 0.0 < R < 5.0 kpc as dashed orange (mass outflow rate) and dotted orange (energy outflow rate) lines to the top left panel.

Figure 14 .
Figure14.Mass loading (top) and energy loading (bottom) as a function of star formation rate surface density split into hot gas with T > 5 × 10 5 K (red dots) and warm gas with T < 5 × 10 5 K (blue dots) at two heights at 0.5 kpc (left) and 1 kpc (right).The red and blue lines represent a best fit in log-log space.