Quokka-based Understanding of Outflows (QED) I. Metal loading, phase structure, and convergence testing for Solar neighbourhood conditions

Multiphase galactic outflows, generated by supernova feedback, are likely to be more metal-rich than the interstellar media from which they are driven due to incomplete mixing between supernova ejecta and the ambient ISM. This enrichment is important for shaping galactic metallicities and metallicity gradients, but measuring it quantitatively from simulations requires resolution high enough to resolve mass, momentum and energy exchanges between the different phases of the outflows. In this context, we present simulations of outflows, driven by SN feedback, conducted using \textsc{Quokka}, a new GPU-optimised AMR radiation-hydrodynamics code. This code allows us to reach combinations of resolution, simulation volume, and simulation duration larger than those that have previously been possible, and to resolve all gas phases from cold neutral medium, $T \sim 100$ K, to hot ionised gas, $T \gtrsim 10^7$ K. In this, a first of a series of papers exploring generation and evolution of multiphase outflows from a wide range of galactic environments and star formation rates, we quantify the extent of selective metal loading in Solar neighbourhood-like environments. We explain the selective metal loading we find as a result of the transport of metals within and between phases, a phenomenon we can study owing to the parsec-scale resolution that our simulations achieve. We also quantify the sensitivity of metal loading studies to numerical resolution, and present convergence criteria for future studies.


INTRODUCTION
Almost all elements heavier than helium, i.e. metals, are manufactured in stars.From their point of origin, metals travel far and wide to populate not only the interstellar medium (ISM) within and the circumgalatic medium (CGM) around galaxies, but also the intergalactic medium (IGM) pervading the cosmos.Because the distribution of metals is driven by several physical processes, such as outflows from stellar feedback, inflows of metal-poor gas from the IGM, and mixing of different gas phases, metals act as tracers of these phenomenon.By following the life-cycle of metals as they traverse a galaxy, we can hope to understand the complex process of galaxy evolution.
An important observational result that links metal abundances to large-scale galaxy evolution is the mass-metallicity relation (MZR; Tremonti et al. 2004), a relatively tight correlation between gas phase metallicities and galaxy stellar masses that extends over several orders ★ E-mail:aditi.vĳayan@anu.edu.au of magnitude in stellar mass.The metallicity, expressed in terms of the oxygen abundance of the gas phase, increases sharply with stellar mass at low galactic masses (albeit with a great deal of scatter) and flattens out at larger masses.In dwarf galaxies one possible explanation for this correlation is the preferential removal of metals in supernova-driven outflows (e.g.Peeples & Shankar 2011;Zahid et al. 2014;Christensen et al. 2018;Forbes et al. 2019).Peeples & Shankar (2011) quantify the selective enrichment of galactic outflows relative to the mean metallicity of the ISM in terms of the "metal expulsion efficiency" or "metal loading factor", , which mathematically is ratio of the metallicity of the outflowing gas to that of the gas present in the ISM.Models with  ≫ 1 have the advantage that they do not require the extreme mass loading factors (often ≳ 100) that are required to explain the mass-metallicity relation in dwarfs in more conventional models where  = 1, i.e., where outflow and ISM metallicities are assumed to be the same (e.g.Finlator & Davé 2008;Davé et al. 2012;Lilly et al. 2013).Moreover, radially-resolved models with  = 1 tend to produce relatively steep metallicity gradients in dwarf galaxies, in contradiction to the observed flatness of dwarf gradients (the mass metallicity gradient relation, MZGR; e.g., Belfiore et al. 2017;Mingozzi et al. 2020;Poetrodjojo et al. 2021); by contrast, models including selective metal loading naturally reproduce the MZGR in dwarfs (Sharda et al. 2021a(Sharda et al. ,b, 2023)).
While these analyses hint at the importance of selective metal loading for galactic properties, both direct observational and numerical explorations have been limited.With regard to the former, though observations of outflows in star-forming galaxies are ubiquitous (e.g.Veilleux et al. 2020, and references therein), to date there are only limited observational constraints on the metal-loading (or mass-loading) factor, and then only for a limited range of outflowing gas phase.For the warm ionised phase, Chisholm et al. (2018) find for a sample of seven galaxies over a wide range in stellar mass that the outflow metallicity can exceed the ISM metallicity by as much as a factor of 50 in dwarf galaxies, though this number might be smaller for larger galaxies.They also find that the bulk of the metals in the outflows arise from entrained ISM, rather than direct SN ejecta.However, metal loading factors derived from absorption studies are fraught with uncertainties related to parameters such as geometry and ionization correction.Cameron et al. (2021) circumvent uncertainties related to ionization by directly estimating the electron temperature from auroral lines in the winds of Mrk 1486, an edge-on disc galaxy with biconical outflows (Duval et al. 2016).They find that outflows along the minor axis of the galaxy are enriched with metals compared to the ISM of the disc.It is worthwhile to note that though both Cameron et al. (2021) and Chisholm et al. (2018) conclude that the outflows in Mrk 1486 are metal enriched, their estimates of the degree of enrichment differ by a factor of two, likely indicative of the level of observational uncertainty.
There is also significant evidence for selective metal enrichment in studies of the faster, hotter components of outflowing gas.Martin et al. (2002) and Stevens et al. (2003) find that X-ray emitting gas in the winds of the dwarf starbursts NGC 1569 and M82, respectively, has a higher  to iron ratio than those galaxies' interstellar media, strongly suggesting the presence of incompletely-mixed type II supernova ejecta in the wind.Recently, Lopez et al. (2020) analysed detailed Chandra spectra of M82 to estimate the metal content of the different temperature phases of its multiphase outflowing gas.Consistent with the conclusions drawn from the aforementioned UV studies of the warm ionised phase, and with the earlier X-ray work, they find that the outflows show gradients in metallicity, which differ by temperature.The relatively cooler warm-hot (≳ 10 6 K) phase, traced by O and Ne lines, retains near-Solar abundance similar to that in the disc, while the hot (> 10 7 K) phase, traced by Si and S, is enriched relative to the disc by up to a factor of ∼ 3.5, although these results are at least somewhat sensitive to the fitting procedure (Ranalli et al. 2008;Konami et al. 2011).
Differential metal loading of the phases of the outflowing gas is an idea supported by theoretical works as well.Melioli et al. (2013) study 3D simulation of galactic winds in a starburst system.They follow chemical evolution of the outflows and find that the metallicity of the lower density (higher temperature) phase may be nearly 4.5 times larger than that of the higher density (lower temperature) phase.A recent study by Emerick et al. (2018) that tracks detailed chemical evolution of ejecta from not only type II SN but also AGB stars in an isolated dwarf galaxy also conclusively shows that metal enrichment preferentially takes place in the hotter gas phase.Andersson et al. (2023) and Rey et al. (2023) reach similar conclusions in simulations reaching higher resolutions -the former reaches a maximum resolution of 1.5 pc, but only in the densest regions of the galactic midplane, while the latter has a peak resolution of 18 pc but maintains this resolution well into the outflow.These findings are in line with the prediction from the seminal paper by Mac Low & Ferrara (1999) that smaller galaxies lose nearly all, if not all, the metals from SN ejection.Long-term simulation of such outflows, explored by Melioli et al. (2015) (see also Fragile et al. 2004;Rodríguez-González et al. 2011), further underscore that the metal-enriched SN ejecta escapes into IGM and may lead to its enrichment (Aguirre et al. 2001).
While isolated galaxy simulations offer deep insights into the overall budget of metals between a galaxy and its surroundings, they lack the spatial resolution required to accurately follow the exchange of metals amongst the phases of the outflowing gas -a physical process that occurs at pc or even smaller scales (e.g., Gentry et al. 2019).Capturing this exchange requires resolutions at this level not just in the galactic plane, but far into the outflow (Vĳayan et al. 2020), and is required for reliable simulation-observation comparison, since observations are generally sensitive only to particular gas phases.In this respect, tall-box simulations emulating a smaller patch of a galactic disc, but offering ∼pc scale resolution, are better poised to study metal fluxes (Kim & Ostriker 2017;Li et al. 2017).
Even the best-resolved of the tall-box simulations carried out to date, however, have limited resolution or limited volume.Creasey et al. (2015) study metal enrichment of galactic winds in an aggregate sense through a suite of simulations that explore the parameter space of gas surface density and gas fraction.Their results indicate that the metal loading of winds emanating from a larger galaxy is weaker even though the absolute metallicity of the winds is higher.However, while they achieve 2 pc resolution, their simulation volume extends to only ±500 pc around the galactic midplane, meaning that they cannot study the phase structure of the outflow except very near the plane.Li et al. (2017) use a suite of tall-box simulations to understand the relationships between mass, energy, and metal loading factors and the underlying gas properties.For their fiducial run, which replicates Solar neighbourhood conditions, they achieve a maximum resolution of 2 pc only within 500 pc of the midplane, and the resolution worsens to 8 pc beyond || > 1 kpc.In order to successfully capture these between different temperature phases, simulations with uniformly high resolution are an imperative.
The simulations of metal loading of outflows offering the highest combination of volume and resolution published to date are those of Kim et al. (2020) and Schneider et al. (2020).The former carry out simulations of metals in galactic winds with a uniform resolution of 256 2 × 1792 with spatial scales of 2 or 4 pc per cell, corresponding to simulation regions that extend to either ±1.8 or ±3.6 kpc around the midplane.The latter use 2048 2 × 4096 cells in a 10 × 10 × 20 kpc domain, achieving 5 pc resolution, but with a 10 4 K temperature floor so the simulations omit neutral material.However, the literature to date lacks a rigorous convergence study identifying convergence criteria and demonstrating converged measurements of metal loading of galactic outflows; the most through study to date is that of Schneider et al. (2020), who find that quantities such the phase structure remain unconverged at their highest resolution of 5 pc.Thus it is unclear if the results from any published simulations are converged.
In this paper we overcome the challenge of high computation costs of high resolution simulations by using the state-of-the-art GPUbased adpative mesh refinement (AMR) code, Quokka, which is significantly faster than CPU-based codes (Wibking & Krumholz 2022).We use this code to examine metal loading of winds in a starforming galaxy using uniformly-high resolution tall-box simulations, reaching combinations of numbers of cells and resolution comparable to or better than the best published to date.We distill the simulation results to express the metal-loading of galactic outflows using two metal-loading factors that track different routes of metal enrichment of outflows, and we demonstrate that our estimates for these quantities are converged in resolution.This paper is the first in a series using Quokka to study the transport of metals in outflows, the Quokkabased Understanding of Outflows Derived from Extensive, Repeated, Accurate, Thorough, Demanding, Expensive, Memory-consuming, Ongoing Numerical Simulations of Transport, Removal, Accretion, Nucleosynthesis, Deposition, and Uplifting of Metals (QUOD ERAT DEMONSTRANDUM, or QED for short).This first paper focuses on quantifying metal loading and phase structure in uniformly-resolved (i.e., non-AMR) simulations of Solar neighbourhood conditions, and on testing for convergence in these simulations.Subsequent papers will explore the use of adaptivity in wind simulations, and will extend the study to other galactic environments.
The remainder of this paper is organised as follows.In Section 2 we describe our numerical methods.In Section 3, we present our primary results regarding both the physics of metal loading and the numerics of achieving converged measurements of it.

Simulation setup and initial conditions
We conduct 3D HD simulations using Quokka (Wibking & Krumholz 2022), an adaptive mesh refinement (AMR) radiationhydrodynamic code optimised for GPUs.For the purposes of this paper we do not include radiative transfer, only radiative cooling.For its hydrodynamic step, Quokka solves the Euler equations of compressible gas dynamics using a method of lines formulation with an RK2 update that is second-order accurate in time and space.We also include gravitational forces provided by a static potential representing the stars and dark matter in a galactic disc; in the present work we do not include gas self-gravity.
All our simulations take place in a 1 × 1 × 8 kpc 3 domain, with the longest dimension along the −axis and the galactic plane centred at  = 0. We use periodic boundary conditions in the  and  directions, and first-order extrapolation boundary conditions in the  direction.For this first study we do not use the AMR capability of Quokka in order to ensure that our simulations have uniformly high resolution throughout the wind and to make testing for convergence straightforward; we will extend this study to AMR simulations in a future work.For this paper, our resolution is uniform throughout the volume, and our cell sizes range from Δ from 32 to 2 pc, corresponding to resolutions from 32 2 × 256 to 512 2 × 4096 cells.We summarise the properties of the simulations in Table 1, and for convenience from this point on we refer to the simulations as FG, where  is the resolution in pc and FG indicates that we are using fixed (i.e., nonadaptive) grids.For comparison, our highest resolution case, FG2, contains ≈ 10× as many cells as the simulations of Kim et al. (2020).It contains a factor of 4 fewer than the simulations of Schneider et al. (2020), but has 2.5× higher resolution and runs for a substantially longer time (≈ 120 Myr as opposed to 70 Myr), allowing it to reach statistical steady-state.
The initial density and pressure profiles are adapted from the Solar Neighbourhood model of TIGRESS simulations (Kim & Ostriker 2017).In this model, the initial density profile follows a double exponential representing a two-phase medium, where  1,2 are the sound speed of the two phases, here set to 7 km s −1 and 70 km s −1 , respectively.The midplane densities are  1,0 = 2.85 H cm −3 and  2,0 = 10 −5  0,1 .The external gravitational potential, Φ ext , is set by the dark matter halo potential and a stellar disc.The dark matter potential is adapted from Kuĳken & Gilmore (1989), and the total potential from the dark matter halo and the stellar disc (reproduced from Kim & Ostriker 2017) is, Here, Σ * = 42 M ⊙ pc −2 ,  * = 245 pc,  dm = 6.4 × 10 −3 M ⊙ pc −3 and  0 is the Galactocentric radius of our simulation box, which we set to be 8 kpc.With these choices of  1,0 ,  2,0 ,  1,2 , and Φ ext , the initial gas surface density, Σ gas , is 13 M ⊙ pc −2 As noted above, while we do not include radiative transfer, we do include radiative heating and cooling.We implement these using a custom cooling source term that is similar, but not identical, to that used by the Grackle library (Smith et al. 2017).We cannot use the Grackle code itself because Grackle does not run on GPUs.Instead, we adopt the tabulated primordial and metal line heating and cooling tables that are included with Grackle, and we re-implement the tabular interpolation routines, as well as the terms for photoelectric heating and Compton cooling (which are not included in the tables themselves) with a temperature floor of 100 K.However, unlike Grackle, we do not include an X-ray heating term.We then integrate the cooling function in each cell in an operator-split manner using Quokka's adaptive RK integrator that runs on GPUs.The model used in Quokka and Grackle includes photoelectric heating at a rate that matches that in the Solar neighbourhood, and produces an atomic medium with well-defined warm and cold phases whose properties are in reasonable agreement with observations; see Appendix A of Wibking & Krumholz (2023) for a detailed discussion of the cooling model and a comparison between it and other commonly-used approaches.

Implementing supernova feedback
We implement SN feedback by injecting thermal energy equivalent to a single SN event, 10 51 erg, into a single cell in the simulation domain.We also add a fixed density of passive scalar Δ SN / cell to the same cell, where  cell is the cell volume, representing metal injection due to SNe; note that the value of Δ SN is arbitrary, as we will discuss in Section 2.3.The total number of such feedback events is determined by the star formation rate corresponding to the initial gas surface density.We use a star formation rate density of 6 × 10 −3 M ⊙ kpc −2 yr −1 which is the SFR used by Li et al. (2017) for Solar neighbourhood conditions.For a Chabrier (2001) IMF, the corresponding surface rate density of SN events is Σ SN = 6 × 10 −5 kpc −2 yr −1 , and the total SN rate in our simulation box is therefore SNe are distributed randomly in the  −  plane and in the −direction their distribution follows a Gaussian with a width of 150 pc.Thus the SN probability per unit volume per unit time is where ℎ = 150 pc and  = Σ SN ℎ −1  −1/2 , so that ∫ P  = Σ SN .For this probability density, the expected number of SNe per time step  in a volume  cell is In practice we implement SN as follows: for each time step of size  on the coarsest AMR level, we compute the number of SNe that will Name Δ (pc)  10); ( 5) & ( 6) Steady state value of the metal loading factor  for  bg = 0 or  O,⊙ (Equation 9).For columns ( 4)-( 6), the central value we report is median over times from 100 − 116 Myr (after steady-state has been established), and the super-and subscripts indicate the temporal 84−th and 16−the percentile, respectively.
occur during that time step by drawing from a Poisson distribution with expectation value Γ SN .For each SN that is to go off during this time step, we determine the position by drawing from a uniform distribution in the  and  directions and from Gaussian distribution with width ℎ in the  direction.We then add thermal energy and passive scalar to the cell enclosing the coordinates for that SN.
The SN events occurring in the initial few Myr of the simulation run create a hot phase.Because this is a volume-filling phase, subsequent SN events are more likely to occur in a low density region of the disc, ensuring that the snowplow radius of the supernova remnant is larger than the resolution of the simulation (Forbes et al. 2016).We emphasise that we do not provide any sub-grid treatment of SN feedback, e.g., injecting radial momentum or kinetic instead of thermal energy.Such models are unnecessary at the resolutions we reach, and are undesirable because they are necessarily resolution-dependent, which can make it impossible to test for or achieve convergence.

Quantifying metal loading
In this section our goal is to explain how we quantify the degree of metal loading in our simulations.

Computing the metal abundance
Before we can compute metal loading, we must first calculate the metallicity in each cell of the simulation.This will be determined by two factors: the initial metal abundance present at the start of the simulation, to which we refer as  bg , and the metal mass added per SN, to which we refer as Δ  ,SN .Note that, thanks to the setup of our simulation, we are free to choose these factors ex post facto -that is, since SNe inject only passive scalar and not mass, and since we know both the total mass density and the passive scalar density in each cell, we are free to compute the metallicity in the cell as an arbitrary linear combination of these two densities after the simulation has already run.Specifically, we write the metal density in every cell as where  is the total mass density and   is the passive scalar density.Note when a SN goes off in a cell the scalar density   in that cell increases by Δ SN / cell , so   increases by Δ  ,SN / cell , as it should.

SN metal injection rate Γ SN ΔM Z,SN
Total outflow metal flux A schematic depicting the meaning of our two metal loading factors,  (Equation 10) and  (Equation 9), and their relationship to the metal injection rate and outflow metal flux.
In practice we choose values of  bg and Δ  ,SN appropriate for oxygen.The oxygen output of a single type II SN is ≈ 1 M ⊙ (Nomoto et al. 2013), and we therefore adopt Δ  ,SN = 1 M ⊙ .We consider a range of values of  bg , which we parameterise in terms of the Solar oxygen abundance:  bg = (/ ⊙ ) O,⊙ , where  O,⊙ = 8.6 × 10 −3 (Asplund et al. 2009).Below we consider values of / ⊙ from 0 to 2.

Definitions of metal loading factors
Now that we have defined the metallicity   in each cell, we are in a position to define the key quantity that we wish to extract from our simulations, the "metal loading factor", generally taken to be the ratio between the metal flux and the star formation rate in the galaxy (see for example Equation 7 of Li et al. (2017), and also Kim et al. (2020), for a definition based on fluxes).However, such a definition assumes that all the metal outflows are solely sourced from the SN activity.A galaxy goes through multiple episodes of star formation over its lifetime, each of which adds to the ambient metallicity of the ISM gas in which future SN events will occur.These events in turn will entrain some of the previously enriched ISM into outflows.Under such circumstances, the metal outflows comprise contributions from not only direct SN ejecta but also the entrainment of enhancedmetallicity ISM, and we wish to define metal loading factors that can track these two channels independently.
Figure 1 provides a schematic illustration of the picture that motivates our definitions: SNe inject metals at a rate Γ SN Δ Z,SN into the region close to the midplane of the galaxy which lies at the bottom of the figure.These metals can end up in either of two boxes labeled "ISM" or "Outflow".Metals that end up in the ISM box are mostly retained by the galaxy and enhance its overall abundance.However, as gaseous outflows are established in the galaxy, some of these retained metals may be entrained with the mass outflows, contributing an amount ⟨⟩  to the total metal flux, where ⟨⟩ is the mean metallicity of the ISM gas being entrained, and  is the mass outflow rate for entrained ISM gas.By contrast, metals produced by SNe that are ejected directly into outflows do not mix with ISM and instead escape the disc -this process is represented by the "Outflows" box.The net metal outflow rate   includes both the entrained ISM and direct outflow components.By contrast, the total mass flux is overwhelmingly dominated by the entrained component, since direct SN ejecta carry very little mass; i.e., the total mass outflow rate is simply  out , the same as the outflow rate for entrained ISM.
To account for the different contributions to the metal outflows, we introduce two different factors to quantify metal loading, viz,  and , which we define as follows.First, we define the net metal outflow rate through a plane of fixed height  at time  as where   is the −velocity of the gas and   is the metal density.
We note here that the net outflow rate as given here includes both inflowing and outflowing gas.However, once steady-state has been achieved most of the gas is outflowing.As a result, the results do not change substantially if we instead consider only outgoing gas.
Similarly, the mass outflow rate through the surface at height  is where  is the total mass density.Since our system is symmetric about  = 0, at least statistically, in practice we will always use [   (, ) +   (−, )] in place of simply   (, ), and similarly for , i.e., we should always understand that when we write   or  at a given  and , the quantity we intend is actually the sum of the fluxes through the + and − surfaces.
We then separately estimate the contribution to the metal flux from entrained ISM as ⟨⟩ , where ⟨⟩ is the average metallicity of material bounded between − and +, i.e., Using these equations, we can define our metal loading factor  as the ratio of the metal outflow rate to the outflow rate that would be expected if the outflows consisted purely of entrained ISM, i.e., For example,  = 2 corresponds to a situation where the metal flux in the outflow is twice what would be expected if the outflow consisted purely of entrained ISM.This quantity is analogous to the  factor defined by Peeples & Shankar (2011).By contrast, the factor  quantifies the fraction of SN metal output that is directly added to outflow without ever mixing with the ISM.We define this quantity as Here the numerator can be interpreted as the metal outflow rate subtracting off the contribution from entrained ISM, while the denominator is the total metal injection rate by SNe.Thus for example a factor  = 0.3 corresponds to a case where 30% of the metals injected by SNe are never mixed into the ISM, and are instead lost promptly; this quantity is analogous to the SN yield reduction factor introduced by Sharda et al. (2021a).

Qualitative simulation outcomes
We begin by describing the qualitative outcome of our simulations in order to orient the reader for the quantitative analysis that follows.
For this purpose we make use of run FG2 evaluated with  bg = 0, though we note that the qualitative behaviour is the same in all runs, and that for phenomena where the value of  bg matters we will show multiple sample values.As SN feedback begins in the system, hot bubbles develop around the injection sites.Within a few Myr, the individual bubbles expand and break out of the disc.SN feedback produces a volume-filling hot gas and subsequent SN explode into this medium.Disc-wide outflows are set up in the galaxy and the initially stratified medium turns multiphase.
Figure 2 shows slices of gas density (), metal density (  ), metallicity (  /), and temperature and the column density along the  axis at a time when steady outflows have been set up in the galaxy.In the outflowing gas, we identify the warm, dense gas likely lifted from the disc.The hotter parts of the outflows, arising from direct injection of SN, are metal-enriched while the cooler parts are comparatively metal-poor.As these different phases propagate out of the disc, they mix and produce regions of metal-poor warm gas surrounded by hot metal-enriched gas; these features are seen particularly clearly in the inset panels, which zoom in on some of the cool clouds.From the column density we can make out the disc comprising cool, dense gas.
We plot the mass and metal outflow rates, as computed from Equation 6 and Equation 7, through the || = 1, 3 kpc surfaces as a function of time in Figure 3.To reduce noise the outflow rates are averaged over a thickness of 5 cells both above and below the surfaces.Both mass and metal fluxes rise initially as individual superbubbles break out and outflows escape from the disc, with the rise occurring first at || = 1 kpc and then later at || = 3 kpc.After ∼ 100 Myr of evolution sustained outflows of mass and metals are set up in the entirety of the simulation domain.The system achieves a near steady-state around this time as subsequent outflow rates fluctuate only at a factor of ≲ 2 level, although there is a slow secular decrease in the mass outflow rate due to the loss of gas mass from the simulation domain through the boundaries at  = ±4 kpc.The mass and metal fluxes through the 1 kpc and 3 kpc surfaces are very similar, as expected given that we are plotting net mass fluxes.However, even if we plot outward-going only mass fluxes, the results are not substantially different, indicating that most of the material that reaches a height of 1 kpc also reaches 3 kpc.The total steady-state mass outflow rate,  ≈ 0.5 − 1 × 10 −2 M ⊙ yr −1 , is comparable to the star formation rate that corresponds to our chosen supernova rate,  * = 6 × 10 −2 M ⊙ yr −1 .Thus the overall mass loading factor in our simulation is ≈ 1.

Bulk loading factors
Once a steady-state has been established, we use the time-averaged outflow properties to estimate the bulk metal loading factors - and  -using Equations 9 and 10, respectively.We compute both quantities as a function of height, and all the values we discuss are averaged over the time interval 100−116 Myr, after the outflow mass and metal fluxes have settled to steady-state.We report these values for every run in Table 1.We note here that we use "net" outflow fluxes, that include both outflowing and inflowing material, but that the results do not change substantially if we use outflowing material only.

Metal loading as a function of height and background metallicity
measures the relative enrichment between the outflowing gas and the gas which has been entrained from the ISM in the outflows.In Figure 4, we show  and  as functions of height from the mid-plane for FG2.We use net outflow rates averaged over 20 pc slabs around each value of .For the case  bg = 0, we find  ≫ 1 indicating that metal outflows are dominated by the highly metal-enriched SN ejecta. decreases at large distances from the mid-plane mostly because of increasing metallicity of the entrained gas, since outflow rates do not change significantly with height (see Figure 3).In contrast to the metal loading factor , the yield reduction factor, , quantifies the proportion of metals added by the SN feedback that are immediately lost to outflows.That  remains high even at large heights (barring the decrease towards the edge of the box which we   9) and  (Equation 10) for a pristine background,  bg = 0.  ≫ 1 indicates that metals outflows are dominated by fresh SN ejecta rather than entrained ISM gas, while  values close to unity indicate that most of the metals added to the galaxy by feedback are lost to outflows.
believe is a result of stochasticity in the simulations) suggests that most of the SN ejected metals might escape the disc and eventually contaminate the CGM.
Because the metal outflow rate and the average metallicity also depend on the level of background enrichment, we expect  to change with  bg .In Figure 5 we show the variation of both  and  with  bg at two different heights.As expected,  decreases with increasing background metallicity because as  bg increases entrained ISM contributes an increasingly large fraction of the outflowing metal flux.At  bg =  O,⊙ , we find  ≈ 1.3, which corresponds to the metal outflows containing a slightly sub-dominant contribution from direct, unmixed SN ejecta and a stronger contribution from metals entrained from the background ISM.Examining the dependence of  on  bg more broadly, we find that, for outflows typical of Solar Neighbourhood conditions, entrained ISM and direct SN ejecta contribute approximately equally for a background metallicity  bg ≈ 0.5 O,⊙ with direct ejecta dominating at lower metallicity and entrained metals at higher metallicity.
We expect that the yield factor should not depend on the ISM  The values of all quantities are normalised to the results of run FG2, the highest-resolution run.Note that  is converged even at 32 pc resolution in the case of an enriched background, but does not converge until much higher resolution for the case  bg = 0.
metallicity  bg , and Figure 5 confirms this expectation:  is nearly independent of  bg , and, as Figure 4, is nearly constant with height as well.A critical conclusion to draw from Figure 5 is that  is quite close to unity, meaning that a significant majority of SN-injected metals are lost to outflows rather than mixing with the ISM.

Testing Convergence
Thus far we have focused on results from run FG2, our highest resolution run.However, we have not yet established that our results are converged at this resolution, and we have yet to establish convergence criteria.We do so by using the metal loading factor ; since this in turn depends on the total mass outflow rate, the metal outflow rate, and the mean metallicity, this implies that we aim for convergence in all these quantities.We conduct a series of runs by progressively increasing the base resolution of the grid and establishing steady state in the outflow rates.We expect to reach convergence eventually because we implement feedback as pure thermal energy without a subgrid recipe.Because the quantities related to metal distribution, i.e., ⟨⟩ and   , depend on mixing between hot and warm phases, once the interfaces between two phases are resolved, we should achieve convergence.
Figure 6 shows the variation of the metal loading factor and the yield reduction factor with resolution; we also report the numerical values shown in the figure in Table 1.Apart from these factors, we also show how the mass and metal outflow rates and the average metallicity change with resolution.All quantities shown are temporal averages of the spatial averages described in Section 2.3 across a 20 pc slab around the height of 1 kpc, though other heights yield qualitatively similar results.Because   and ⟨⟩ depend on background metallicity, we show the variation of these quantities for both pristine (left) and Solar enriched (right) backgrounds.For  bg = 0, between the lowest and highest resolutions, the mass (metal) outflow rates increase (decrease) by a factor of ≲ 2 as we go from 32 pc to 2 pc resolution (  = 256 to 4096 cells in the  direction).In the same interval, the average metallicity suffers a steeper decline, resulting in a factor of ∼ 5 increase in .This points to the importance of resolving the interfaces between the temperature phases where metal exchange primarily occurs -the mean ISM metallicity ⟨⟩ is lower in the  bg = 0 runs at higher resolution because increasing resolution leads to less numerical mixing between the hot and cold phases, and thus to less metal contamination of the cold gas.The left panel of Figure 6 shows that though convergence may been achieved in the mass outflow rate at relatively modest resolution, the metal loading factor may not necessarily be converged as a result of this effect.Consequently,  does not appear to converge until ≈ 4 pc resolution.By contrast the yield reduction factor, , depends only on the metal outflow rate for  bg = 0 for which  ≫ 1.Therefore, its convergence curve follows that of the metal outflow rate.
A high background metallicity erases almost all variation in ,   , and ⟨⟩ with resolution, such that it appears that these quantities are converged even at 32 pc resolution.There is but slight variation in   of factor ≲ 0.5 which translates into similar variation in .We stress here that for an enriched background, it is easier to achieve convergence in metal outflow rates and consequently the metal loading and the yield reduction factors, simply because at higher background metallicity numerical diffusion from the hot phase into the cool ISM represents a smaller perturbation.
Examining Figure 7, we see that all these phases are populated both close to the disc and far from it, but that the relative contributions vary with height.Both gas and metal masses in the region closer to the disc are dominated by CNM and WNM, while the CNM is much sparser in the extraplanar regions.In the || > 1 kpc region, the bulk of the mass lies in the WNM and WIM, though this region also hosts more WHIM and HIM than the region closer to the disc.We note here that within  < |1| kpc HIM amd WHIM host more metals relative to the gas mass they carry.In the extra-planar regions, HIM and WHIM host most of the metals.For HIM this shows that most of the metals do not mix with the ISM and quickly escape into extra-planar region.WHIM acquires metals by cooling of HIM and heating of cooler phases by means of mixing.In the remainder of this section we examine the properties of the outflow as a function of phase; because CNM and UNM are a subdominant (but, we emphasise, not completely negligible) component in the outflow region, for simplicity in the remainder of this section we will group these phases together with WNM as a single neutral phase.

Phase distribution by mass and flux
While Figure 7 shows the distribution by mass, it is interesting to contrast this with the distribution by flux.To explore this, we show the distribution of mass (top) and metal (middle) fluxes passing through 2 pc thick slabs at different heights in Figure 8.For comparison, in the bottom panel of Figure 8 we show the partial density in each phase, defined as the total mass of material in that phase divided by the volume containing all material (as opposed to the mass of each phase divided only by the volume it occupies).
The phase structures of mass and metal fluxes are quite different.Focusing first on the former, the neutral, WIM and WHIM phases carry ∼ 10% of mass flux at 1 kpc but this increases to 50% at 4 kpc; by contrast, that the neutral phase dominates the mass budget at all heights.Given that the cooler phases are responsible for most of the mass at all heights, and that the increase in their partial densities with height is much smaller than the increase in their contribution to the flux, we can conclude that these phases are being accelerated, rather than forming via condensation of hot gas.We note here under ballistic conditions we should expect a decrease in momentum of the cooler phases, neutral in particular, simply due to the inertia it carries.In fact, we find the opposite, a point whose significance we discuss in Section 4.
With regard to metals, the key conclusion to be drawn from Figure 8 is that most of the metal flux is carried in the hot phase of the outflows, though the proportion changes with height.Closer to the midplane, at 1 kpc, hot outflows carry nearly 90% of the metal flux, while at the larger height this decreases but remains > 50% of the total.This is in marked contrast to the distribution of metals by mass shown in Figure 7, where even at || > 1 kpc most of the metal mass lies in gas with  ≲ 10 4 K.This difference can be attributed to the velocity structure of the gas: while the hot phase contains less metal mass, its outward velocity is much larger, and thus it carries a majority of the metal flux.The situation is quite different for the mass flux, which is predominantly hot at || = 1 kpc, but where the balance shifts in favour of warm (T≲ 10 4 K) as gas moves from 1 kpc to 4 kpc.In regions closer to the midplane, cold and dense phase dominates the mass content, while for the outflowing gas mass shifts to higher temperatures.

Phase-wise metallicity
Thus far we have examined the distribution of mass and flux with respect to phase; however, these quantities are not easily accessible via observations.Instead, what observations can probe is differences in the abundance of one phase of the outflowing gas versus another, and differences with respect to height within a single phase.We investigate these differences by computing time-averaged abundances (Equation 8) over 1 kpc-wide regions of the simulation domain at times after 100 Myr, once steady-state has been established.We compute these averages separately for each of the phases listed above, though for convenience and to avoid cluttering our plots we again group all the neutral phases (CNM, UNM, WNM) together for this purpose.
We plot the average metallicities as a function of distance from the mid-plane for each phase in the top row of Figure 9; the three columns show three different background metallicities,  bg / O,⊙ = 0, 0.2 (comparable to the metallicity of the Small Magellanic Cloud and of dwarf starbursts whose outflows have been studied in observations), and 1.We see that the hot phase, carrying the fresh SN ejecta, is the most metal-rich at all heights, but that the difference between its metallicity and the mean metallicity of the cooler phases is a function of both height and background metallicity; close to the mid-plane and for  bg = 0, the hot phase is as much as ≈ 5× more metal-rich than any of the cooler phases, while for  bg =  O,⊙ , the difference drops to at most tens of percent.At an intermediate metallicity of  bg = 0.2 O,⊙ , the difference is a factor of ∼ 2 at low heights, dropping to tens of percent by ≈ 3.5 kpc.In general the fact that we see both HIM metallicity decreasing and cool phase metallicity increasing with height indicates that there must be material exchanged between the phases in both directions.That is, some metal-poor cool gas must be evaporating into the hot phase in order to explain why the mean HIM metallicity decreases, while some hot material must be condensing into the cooler phases to explain why the WIM, WHIM, and neutral metallicity increases.In the case  bg =  O,⊙ where the neutral and hot phases are at similar mean metallicity even at  = 0, this exchange has little effect on either phase, while its effects are much more dramatic for the case  bg = 0, where the hot and neutral phases have very different metallicities near the mid-plane.However, recalling the discussion of mass and metal fluxes in the previous section, we remind readers that, while this exchange occurs, it is not enough to significantly alter which phases carry the bulk of the metal flux.
While the difference in metallicity between phases is of obvious interest from the standpoint of physical interpretation, it is also difficult to probe with observations due to the challenge of cross-calibrating absolute metallicities between, e.g., X-ray and optical data.For this reason, an alternative quantity that is often reported instead is the metallicity of the outflow relative to the ISM, or the variation of metallicity with distance from the galactic plane, within a single phase.The bottom panel of Figure 9 identical to the top but instead of showing the absolute metallicity, we show metallicity normalised to the mid-plane value.Predictably, for pristine backgrounds, the neutral phase experiences the highest enrichment with respect to its mid-plane metallicity, while the hot phase becomes more metal-poor, albeit by a smaller amount.The relative change in the metallicity of the cool phases is ≈ 50% for  bg = 0.2 O,⊙ , and this drops to ≲ 10% for  bg =  O,⊙ .The latter is likely unmeasurable at the accuracy of current metallicity diagnostics, but the former may well be observable, and indeed may already have been observed, a point to which we return below.

Implications for the origin of the mass-metallicity and mass-metallicity gradient relations
The strongest result from our work is Figure 5 showing that the yield reduction factor  ≈ 0.8 for star formation in Solar neighborhood-like conditions, i.e., ≈ 80% of the metals in SN ejecta are never mixed with the surrounding ISM and instead escape with the outflowing gas.An important subsidiary point is that this is true despite the fact that the metal loading factor  that characterises the ratio of outflow metallicity to ISM metallicity is relatively modest for Solar metallicity backgrounds.That is, the fact that  ≲ 2 for a metal-rich galaxy does not mean that metal loading is modest, as some authors have argued (Kim et al. 2020).Instead, high  coexists with moderate  simply because for a very metal-rich background and moderately mass-loaded winds, even very high SN metal loss only enhances wind metallicity mildly compared to ISM metallicity.Quick expulsion of SN-produced metals is a key piece of physics for understanding the mass-metallicity relation (MZR) and the massmetallicity gradient relation (MZGR), particularly in dwarf galaxies (Sharda et al. 2021a,b).Analytical models for these relations find that the observed relatively flat MZGR at low galaxy masses can be understood only if  is near unity, as we find.Other semi-analytical models and analysis of cosmological simulations have drawn similar conclusions about metal retention in the ISM of dwarfs (Ma et al. 2016;Pandya et al. 2021).Our work strongly supports the view that the paucity of metals in the ISM of dwarf galaxies is a direct result of heavy metal loading of outflows coupled with moderate mass loading (e.g., Forbes et al. 2019), rather than extreme mass loading as has been suggested elsewhere (Davé et al. 2011).
Preferential ejection of SN-produced metals has also been invoked to explain how the CGM of star-forming galaxies came to hold nearly as much oxygen as the disc (Tumlinson et al. 2011;Peeples et al. 2014).Tumlinson et al. (2011) posits that in order to enrich the CGM to the observed levels over a reasonable time scale, most, if not all, of the oxygen produced in SN should be carried out by the outflows, though much of this material is re-accreted over ∼ Gyr timescales and participates in future star formation episodes.This picture is consistent with the narrative set by Figure 8 that metal flux balance continuously shifts towards warmer phases which may not possess the momentum required to eventually escape the galaxy's potential.

Implications for the phase structure of outflows
A second important conclusion to draw from our work is with regard to the phase structure of outflows.The outflow we produce in our simulations has a structure similar to that proposed by Thompson & Krumholz (2016) and Krumholz et al. (2017), whereby the cooler phases exist throughout the outflow and are primarily the result of acceleration of pre-existing cool gas out of the plane, rather than re-condensation of hot gas off the plane (e.g.Thompson et al. 2016;Schneider et al. 2018).The survival of these cool clouds in the outflow in our simulations appears to be a result of efficient radiative cooling,  8) in 1 kpc-wide sections of the simulation domain, separated by temperature phases.The top row shows metallicity on an absolute scale, while the bottom shows metallicity normalised to the mean metallicity of each phase at  = 0. similar to the effect seen in other simulations with similarly-high spatial resolution (e.g., Schneider et al. 2020;Kim et al. 2020).Such a picture is consistent with recent high spatial resolution observations of outflows, which show that the neutral and molecular phases are present even at the outflow base (e.g.Leroy et al. 2015;Martini et al. 2018;Noon et al. 2023), and kinematics suggesting that neutral material accelerates with distance from the galaxy (Yuan et al. 2023).As a result of this structure, neutral or cool ionised material dominates the outflow mass at all heights, but because it accelerates slowly it only becomes a major constituent of the mass flux at heights ≳ 2 kpc.
Metals add an important dimension to this picture because, especially in the case  bg = 0 where the metals are initially present only in the hottest phase, they effectively act as Lagrangian tracers of exchange between phases.The story told by these Lagrangian tracers is that exchange between the HIM phase and the cooler phases is not zero, but is surprisingly small.In particular, recall from Figure 8 that at no height does the HIM contribute more than ∼ 10% of the mass.Thus even if only ∼ 10% of the essentially metal-free cooler gas were to evaporate into the hot phase, this would be sufficient to dilute its metallicity down by a factor of two.In fact, Figure 9 shows that the decrease in mean metallicity of the hot phase is considerably smaller than this, meaning that only a few percent of the initially-cooler material can be added to the hot phase over the 4 kpc distance that we track the outflow.Conversely, the fact that the metal flux in the hot phase decreases by a factor ≲ 2 between 0.5 and 4 kpc (cf. Figure 8) implies relatively small loss of hot, metal-rich gas into the cooler phases.Thus the basic picture toward which we are driven is one where, at least out to 4 kpc, the different gas phases for the most part maintain their identities.There is substantial exchange of momentum, as is required to accelerate the cool gas, but not a great deal of exchange of material.
Of course we emphasise that these conclusions apply only up to 4 kpc.A number of simulations with lower resolution but larger volume (e.g., Schneider et al. 2020), as well as analytic models (e.g., Fielding & Bryan 2022), suggest that there should be more exchange between phases at larger heights.We cannot rule out this possibility.However, we also caution that these conclusions are based on simulations and models that do not allow cooling past 10 4 K, and thus the cool material that becomes hot in these models is assumed to all be in the form of WIM.In fact, we find that even at a height of 4 kpc neutral material represents an equal contribution to the mass flux, and a dominant contribution to the total mass; as noted above, observations of M82 support this conclusion.It is therefore unclear to what extent these models are applicable.

Comparison with other theoretical works
Qualitatively our results reinforce some of the general conclusions of earlier theoretical works, although the physics implemented in each of these is somewhat different, and though our exploration of the effects of varying the background metallicity allows us to draw somewhat different conclusions.Our simulation setup most closely resembles those of Li et al. (2017) and Kim et al. (2020), and our SN injection recipe is similar to Li et al.'s as well, albeit without the contribution of SN Type Ia.Li et al.'s metal loading factor reduces to our definition of  in the case of metal poor backgrounds (see Equation 10 in case of  ≫ 1) and we find that both quantities are similar in value -although it should be noted that the box size used in Li et al. (2017) is smaller than ours in volume.
Compared to Kim et al. (2020), we lack their self-consistent treatment of star formation, but we reach higher resolution over a larger volume.Overall our estimates of the metal outflow rate (Figure 3) and the enrichment ratio, ⟨⟩/⟨ ISM ⟩ (Figure 9), are in qualitative agreement with their "R8" model.In particular, we find that the hot phase carries most metals, that the metal loading factor  is modest ( ≲ 2), and an analogous plot for Figure 8 for a Solar metallicity background shows that at increasingly higher altitudes the metal flux resembles the mass flux as the metal contribution from the entrained gas increases.However, Kim et al. (2020) do not consider the fraction of SN-injected metals lost to outflows, , which is arguably more important from the standpoint of the MZR and MZGR than the metal loading factor .They therefore do not reach the two critical conclusions we reach, namely that most SN-injected metals are lost, and that this loss implies much higher metal loading factors in galaxies with lower mean ISM metallicities.
In addition to previous tall box simulations, we can also compare our results to previous simulations of isolated dwarf galaxies.For such a setup, Emerick et al. (2018) follow detailed chemical evolution of several different ion species and evolve for long timescales, albeit a much lower resolution than we achieve.Their results also point towards poor metal retention, with nearly ∼ 90% of the metals generated by SN feedback being lost, though some of these lost metals may be re-accreted in future.Schneider et al. (2020) also simulate an isolated dwarf, though their simulation follows a starburst galaxy with a much higher star formation rate per unit area than our Solar neighbourhood conditions.Interestingly, they appear to find considerably more rapid phasemixing than we do.They report that the concentration of the passive scalar that they inject into their hot phase is diluted by more than a factor of two even within 2 kpc of the galactic plane, whereas we find lesser dilution even out to 4 kpc (c.f. Figure 9).The cause of the difference is unclear.One obvious candidate is resolution, since at their resolution of 5 pc we find that  is still not fully converged, and lower resolution promotes mixing.However, there are other possible explanations as well, including the differences in star formation rate, initial conditions, and problem geometry between the two simulations, and differences in the hydrodynamic scheme, to which mixing can be sensitive -in particular, their scheme uses PLM reconstruction, which is lower-order than the PPM method we use in Quokka, and thus is likely to produce stronger numerical mixing.

Comparison with observations
Results from our work can be directly compared with observations in both optical/UV and X-ray bands.For the former, Chisholm et al. (2018) measures the metallicity of outflowing gas in galaxies covering several orders of magnitude in mass.Since they use UV absorption, they are able to trace the phases closest to the WIM and WHIM phases described in Section 3.3.2.Their measurements support the conclusion that the outflowing gas is heavily metal loaded with respect to the host's ISM, with the amount of metal loading being larger for metal-poor dwarf galaxies than for more metal-rich galaxies.This is at least qualitatively consistent with our findings, in particular Figure 9, where we find that the for WIM and WHIM the metallicity is larger in outflowing gas (|| ≳ 1 pc) than in midplane gas, but that the difference decreases as the overall galaxy metallicity increases.Interestingly, their entrainment fraction, by which they estimate the fraction of metals in the outflows arising from entrained ISM gas, is ≳ 0.8.This at might first seem at odds with our conclusions that, at least in dwarfs,  ≫ 2, i.e., direct SN ejecta dominate the outflow; Chisholm et al.'s result corresponds to  < 2. However, the contradiction is resolved if we recall from Figure 8 that WIM and WHIM together carry roughly half the mass flux by only ≈ 10% of the total metal flux, precisely because most of the direct SN-ejected metals are carried in the hot phase do not mix into the WHIM or WIM.Consequently, our simulation is consistent with Chisholm et al.'s conclusion that entrained metals dominate, provided that we recognise that this conclusion is limited to the phases that are accessible via UV spectroscopy, and is not true of the outflow as a whole.This finding thus highlights the importance of combining observations that probe more than one phase.
With regard to X-rays, Lopez et al. (2020) analyse Chandra observations of the outflowing gas in M82.They follow the warm-hot and the hot phases up to a distance of ∼ 3 kpc from the disc, traced by O (along with Ne, Mg, Fe) and Si (and also S), respectively.The abundances of both these phases is nearly flat at 1 − 1.5 times the Solar value outside the central injection radius, identified as 500 pc.Such a trend similar to the  bg =  O,⊙ panel in Figure 9, both qualitatively and quantitatively, though we caution about putting too much weight on this agreement given that our simulation conditions are intended to represent the Solar neighbourhood, not a starburst such as M82.The X-ray surface brightness maps show a steady decline in X-ray luminosity towards regions of higher altitude.Though we plan to compute the X-ray emission properties of our simulations in detail in a later paper, we can predict a similar trend on the basis of Figure 2, noting that gas is hotter closer to the disk.Lopez et al. (2020) fit models to the spectra from different regions of the wind extract temperature of the emitting gas.They find that the spectra in most regions outside the disc are consistent with the presence of gas at two distinct temperatures, ∼ 0.4 − 0.6 keV and ∼ 0.8 − 1.7 keV.As can be seen from Figure 7, we also predict significant amount of extra-planar gas in the HIM, which can be as hot as ∼ 10 7 K (= 1 keV).

Caveats
Some pieces of physics that may affect outflow properties and phase structures are not yet implemented in QED, and we therefore pause here to note these caveats.One is self-consistent star-formation and pre-SN feedback.For instance, it has been suggested that clustering of SNe may alter the mass loading of the outflows (Smith et al. 2021), and that this in turn is regulated by pre-SN feedback (e.g., Jeffreson et al. 2021).Any increase in mass loading will affect the metal loading factors of galaxies with a non-zero background metallicity, and it is conceivable that clustering also directly alters phase mixing.However, we note that this effect is less important in simulations that include pre-SN feedback than in earlier ones omitting it.
A second missing piece of physics is cosmic rays (CRs), which some authors find are capable of driving outflows with factor unity of mass loading (Girichidis et al. 2016;Pakmor et al. 2016;Simpson et al. 2016).CRs drive cooler and denser outflows implying they affect mostly the cooler phases of the outflows (Girichidis et al. 2018), which will be responsible for metal outflows only in non-zero background metallicity.The true importance of CRs for outflow driving is, however, extremely uncertain due to its dependence on poorly-constrained parameters of CR transport (e.g., Crocker et al. 2021a,b).
Finally, all gas in QED cools at a rate identical to gas at Solar metallicity -that is, our cooling is not computed self-consistently with the spatially-varying metallicity.Gas entrained from a sub-solar metallicity ISM should cool more slowly, which might affect the overall phase structure of the outflows.However, it should be noted that in the case of a sub-solar metallicity background, the contribution from this phase to the overall metal loading is also reduced.In future work, we intend to remove the inconsistency as we explore how  and  vary with environment.Upcoming iterations of QED will also include a chemistry network that to follow the evolution of individual species rather than lumping all metals together.

SUMMARY AND CONCLUSIONS
We present results from 3D high-resolution simulations of a patch of a Milky Way-mass galaxy using AMR-based code Quokka, optimised for GPUs.Thanks to this optimisation, QED simulations are able to reach a combination of resolution, volume, and run duration that exceeds any published to date.Our initial setup comprises an initially uniform gas disc with properties modelled on the Solar neighbourhood.Supernova feedback is injected by adding pure thermal energy to cells at a rate consistent with the expect SN rate for the gas surface density, and at random locations drawn from a Gaussian height distribution.We tag the SN ejecta with a passive scalar, representing SN-produced metals, which is then advected with gas.This enables us to track the eventual fate of the metals injected into the ISM by SNe.
We find that the simulation quickly develops large scale outflows, which settle into approximately steady state mass and metal outflow rates after ≈ 100 Myr.We quantify the metal loading of the outflow in this steady state in terms of two dimensionless factors,  and .The former is the classical "metal loading" factor and describes the enhancement in outflow metal flux compared to what would be expected if the outflow consisted purely of entrained ISM (i.e., with no contribution from unmixed SN ejecta), while the latter quantifies the fraction of metals ejected by SN that end up in metal outflows without ever mixing with the ISM.
Our main findings are: (i) The metal loading factor  is greater than unity, meaning that the outflowing gas carries more metal flux than would be expected if it consisted solely of entrained ISM.However, the amount by which  exceed unity depends on the background metallicity of the galaxy.Entrained metals dominate the metal flux ( < 2) for galaxies with metallicities ≳ 20% of Solar, while direct SN ejecta dominate in more metal-poor systems.
(ii) By contrast the yield reduction factor, , which characterises the fraction of SN-injected metals that are lost promptly to the outflow, does not depend on the background metallicity and is fairly close to unity,  ≈ 0.8 − 0.9.Thus most of the metals produced by stars leave the disc of the galaxy promptly.Theoretical models for the origin of the mass-metallicity and mass-metallicity gradient relations favour such large values of , and our simulation results provide physical backing to the large values indirectly inferred from these models.
(iii) The phase distributions of mass, mass flux, and metal flux are all very different in metal-poor galaxies.In such galaxies, metals are mostly carried in the fast-moving hot phase, but this phase supplies only ≈ 50% of the mass flux and ≲ 10% of the mass; instead, most mass resides in neutral gas.These differences are a result of there being relatively little exchange of mass between the phases.However, as the background metallicity increases, the differences between the phase balance of the mass and metal flux diminishes, because as entrained ISM becomes more metal-rich it constitutes a larger and larger share of the total metal flux.
(iv) Variations in the outflow metallicity between different phases at fixed height, and within a single gas phase as a function of height, provide a direct and powerful diagnostic of outflow physics, one for which our simulations make definite predictions.In particular, we find that the hot phase should be more metal-rich than warm ionised or neutral gas at fixed height, and that hot gas metallicity should mildly decrease with height, while warm ionised and neutral gas metallicity increases.All of these effects are magnified at low background metallicity and suppressed at high background metallicity.
(v) Capturing metal loading and the balance of metals between different ISM phases in numerical simulations requires very high resolution.We find that these quantities only converge at resolutions of ≈ 2 − 4 pc, while lower-resolution simulations tend to underestimate the extent of metal-loading.

Figure 2 .
Figure 2. Slices of gas density, metal density, metal abundance, and temperature in run FG2 at time  = 115 Myr from left to right, respectively.The right column shows the column density along the  axis.The horizontal line indicates the initial scale height of the disc.For reasons of space we show only the top half of the simulation domain, but remind readers that the domain extends to −4 kpc below the plane as well.Inset panels zoom in on example regions a high resolution.

Figure 3 .
Figure 3. Mass (top) and metal (bottom) outflow rates through surfaces at  = 1, 3 kpc, computed from Equation 6 and Equation 7 for FG2.The figure shows that, after an initial transient, the outflow rates settle down to nearsteady-state values at times ≳ 100 Myr.

Figure 4 .
Figure 4. Metal loading factors  (Equation9) and  (Equation10) for a pristine background,  bg = 0.  ≫ 1 indicates that metals outflows are dominated by fresh SN ejecta rather than entrained ISM gas, while  values close to unity indicate that most of the metals added to the galaxy by feedback are lost to outflows.

Figure 5 .
Figure 5. Same as Figure 4, except now showing  and  at fixed heights of 1 and 3 kpc but for varying  bg .The horizontal dotted line is at  = 2 indicates where direct SN ejecta and entrained ISM contribute equally to metal outflows; this occurs at  bg / O,⊙ ≈ 0.5.

Figure 6 .
Figure 6.Time-averaged values of metal-loading factors  and , mass outflow rate , metal outflow rate   , and mean metallicity ⟨ ⟩, all computed for a 20 pc slab around |  | = 1 kpc and for  bg = 0 (left) and  bg =  ⊙ (right), as a function of simulation resolution.The values of all quantities are normalised to the results of run FG2, the highest-resolution run.Note that  is converged even at 32 pc resolution in the case of an enriched background, but does not converge until much higher resolution for the case  bg = 0.

Figure 7 .
Figure 7. Mass-weighted, time-averaged temperature-density histograms for gas (left) and metals (right), averaged over the regions |  | < 1 kpc (top) and |  | > 1 kpc (bottom); metal masses are computed for the case  bg = 0.The horizontal lines indicate the temperature thresholds for separating the neutral phase (which combines the CNM, UNM and WNM), the WIM, WHIM, and the hot phase.In regions closer to the midplane, cold and dense phase dominates the mass content, while for the outflowing gas mass shifts to higher temperatures.

Figure 8 .
Figure 8.The phase separated fluxes of mass (top) and metal (middle) and the average mass density (bottom).Though HIM is responsible for most of the mass flux close to the disc, the balance shifts towards the WIM and WNM phases at larger  values.In case of the metal flux, the hot material dominates all the way till the edge of the simulation domain.

Figure 9 .
Figure 9. Mass-weighted metallicities (Equation8) in 1 kpc-wide sections of the simulation domain, separated by temperature phases.The top row shows metallicity on an absolute scale, while the bottom shows metallicity normalised to the mean metallicity of each phase at  = 0.

Table 1 .
Summary of runs.(1) Name of the run; (2) Base resolution of the grid in pc; (3) Number of cells in each direction on the base grid; (4) Steady state value of the metal loading factor  (Equation