Boosting galactic outflows with enhanced resolution

We study how better resolving the cooling length of galactic outflows affect their energetics. We perform radiative-hydrodynamical galaxy formation simulations of an isolated dwarf galaxy ($M_{\star}=10^{8}\, \mathrm{M}_\odot$) with the Ramses-RTZ code, accounting for non-equilibrium cooling and chemistry coupled to radiative transfer. Our simulations reach a spatial resolution of $18 \, \mathrm{pc}$ in the interstellar medium (ISM) using a traditional quasi-Lagrangian scheme. We further implement a new adaptive mesh refinement (AMR) strategy to resolve the local gas cooling length, allowing us to gradually increase the resolution in the stellar-feedback-powered outflows, from $\geq 200 \, \mathrm{pc}$ to $18 \, \mathrm{pc}$. The propagation of outflows into the inner circumgalactic medium (CGM) is significantly modified by this additional resolution, but the ISM, star formation and feedback remain by and large the same. With increasing resolution in the diffuse gas, the hot outflowing phase ($T>8 \times 10^{4} \, \mathrm{K}$) systematically reaches overall higher temperatures and stays hotter for longer as it propagates outwards. This leads to two-fold increases in the time-averaged mass and metal outflow loading factors away from the galaxy ($r=5\, \mathrm{kpc}$), a five-fold increase in the average energy loading factor, and a $\approx$50 per cent increase in the number of sightlines with $N_{\text{OVI}} \geq 10^{13}\, \mathrm{cm}^{-2}$. Such a significant boost to the energetics of outflows without new feedback mechanisms or channels strongly motivates future studies quantifying the efficiency with which better-resolved multiphase outflows regulate galactic star formation in a cosmological context.


INTRODUCTION
Multiphase galactic outflows are ubiquitous across the galaxy population, with their range of molecular, cold, ionized and hot phases detected across the electromagnetic spectrum (see e.g.Veilleux et al. 2020;Laha et al. 2021 for reviews).By removing mass from the dense ISM and heating the gas surrounding galaxies in the CGM, outflows provide a fundamental mechanism to regulate the star formation of galaxies across cosmic time.But a detailed understanding of how galactic outflows launch, how they propagate through the ISM and how they interact with the CGM remains elusive, and a key challenge for modern galaxy formation theories (see Somerville & Davé 2015;Naab & Ostriker 2017 for reviews).
Low-mass dwarf galaxies provide valuable clues to address these questions, as their shallower gravitational potential wells makes them acutely sensitive to their internal stellar processes.In particular, supernovae (SNe) explosions provide an important engine to inject energy and momentum into the ISM and launch powerful galactic winds capable of escaping such low-mass systems (e.g.Chevalier & Clegg 1985;Dekel & Silk 1986;Murray et al. 2005).Despite this well-established picture, large uncertainties however remain in the quantitative efficiency with which these SN-driven outflows regulate star formation in dwarf galaxies.
Already within the ISM, additional stellar processes and feedback channels (e.g.winds, radiation; see Naab & Ostriker 2017 for a review) can modify the properties of the gas in which SNe explode, ★ E-mail: martin.rey@physics.ox.ac.uk leading to qualitative and quantitative changes in the venting behaviour of low-mass galaxies (e.g.Emerick et al. 2020;Agertz et al. 2020;Smith et al. 2021;Fichtner et al. 2022).Moreover, for a given feedback model, the spatial distribution of the energy injection and the clustering in space and time of SNe plays a major role in the ability to build a powerful outflow (e.g.Girichidis et al. 2016;Kim et al. 2017;Fielding et al. 2018;Ohlin et al. 2019), making outflow properties strongly sensitive to the underlying model for where young stars form and spatially distribute (e.g.Andersson et al. 2020Andersson et al. , 2023;;Steinwandel et al. 2023).Further, the amount of feedback energy and momentum available from stellar evolution at the low metallicities of dwarf galaxies remains debated, making the total feedback budget itself an important uncertainty in such systems (e.g.Gutcke et al. 2021;Prgomet et al. 2022).All these factors contribute to the large spectrum of outflow energetics reported in the literature, with the predicted efficiency with which outflows carry mass and energy out of dwarf galaxies spreading over orders of magnitudes at a given stellar mass (see e.g.Muratov et al. 2015;Hu 2019;Emerick et al. 2020;Smith et al. 2021;Pandya et al. 2021;Andersson et al. 2023;Steinwandel et al. 2024 for recent studies).
Observational constraints are thus instrumental to help pinpoint the energetics of galactic outflows in low-mass galaxies, but the multiphase nature of such outflows makes this a challenging task.The different cold, warm and hot gas phases are probed through a diversity of observational techniques, each with their own wavelength window and challenges (see Collins & Read 2022 for a review).Moreover, there is no consensus as to the respective importance of each phase in regulating star formation.The cold-to-warm phase (≈ 10 4 K) is expected to dominate the mass budget of the outflow (e.g.Kim & Ostriker 2018;Kim et al. 2020;Fielding & Bryan 2022), but studies targeting this temperature range, i.e. through H or ultraviolet (UV) absorption lines, have both reported above-unity mass-loading factors (Chisholm et al. 2017;McQuinn et al. 2019;Schroetter et al. 2019) and much lower values (Marasco et al. 2023).Alternatively, the hot outflow phase (≥ 10 5 K) is expected to dominate the energy budget, which could be key to pressurize the surrounding CGM and prevent further star formation by delaying gas inflows onto the galaxy (e.g.Hu 2019;Li & Bryan 2020;Carr et al. 2023).However, this phase is challenging to observe, emitting X-rays that are only accessible for limited samples of low-mass dwarfs (e.g.Heckman et al. 1995;Summers et al. 2003Summers et al. , 2004;;Ott et al. 2005;McQuinn et al. 2018).This makes obtaining a holistic overview of outflow thermodynamics difficult observationally, and robust predictions of the multiphase structure of outflows paramount theoretically.
Beyond the launching mechanics of the wind and its initial structure within the ISM, a key aspect to achieving accurate predictions is to robustly capture the propagation of the multiphase outflow as it escapes the galaxy.A spectrum of processes can indeed efficiently transfer mass and energy between gas phases during the propagation, in particular shocks (e.g.Klein et al. 1994), hydrodynamical instabilities and turbulence (e.g.Gronke & Oh 2018, 2020;Fielding et al. 2020;Kanjilal et al. 2021;Tan et al. 2021), cooling and heating of the different gas phases (e.g.Field 1965;Begelman & McKee 1990;Sharma et al. 2010Sharma et al. , 2012;;McCourt et al. 2012;Voit et al. 2015), the interaction with a surrounding hot or already multiphase CGM medium (e.g.Armillotta et al. 2016;Brüggen & Scannapieco 2016;Gronke et al. 2022), or a combination of the above (see Faucher-Giguere & Oh 2023 for a review).The relative efficiency of each of these processes remains debated as they depend on the local gas and radiation conditions, but they introduce characteristic length scales that should be resolved to accurately capture how the original structure of outflows is reprocessed as they expand away from the galaxy.A key common feature for these processes, however, is that their characteristic length-scales are invariably small (≪ 10 pc), posing a huge numerical challenge to simulations that wish to model entire galaxies (≥ kpc).This issue is further compounded by the Lagrangian behaviour of most galaxy simulations, focussing the computational effort and resolution where gas is dense (i.e. in the ISM) but quickly degrading it in the diffuse gas to speed up the computation.Density and temperature contrasts in the launched outflow close to the galaxy are then numerically smoothed out as it expands away and the resolution worsens.This leads to artificial over-mixing of the gas phases and the suppression of thermal instabilities which we illustrate in Figure 1 (see also the discussion in Hummels et al. 2019).
Figure 1 shows the gas cooling length density-weighted along the line of sight in the edge-on isolated dwarf galaxy used in this study.The cooling length is ≤ 100 pc in the outflowing gas escaping the galaxy (top, blue to white to red), but is marginally or underresolved in a traditional galaxy formation simulation relying on a quasi-Lagrangian resolution strategy (bottom left, blue).The multiphase structure of simulated galactic outflows is thus likely to be under-resolved and numerically suppressed, reminiscent of the numerical limitations in modelling the larger-scale diffuse CGM (e.g.Peeples et al. 2019;Hummels et al. 2019;Suresh et al. 2019;van de Voort et al. 2019).These challenges have motivated subgrid implementations and new hydrodynamical methods to model multiphase galactic winds (e.g.Huang et al. 2020Huang et al. , 2022;;Weinberger & Hernquist 2023;Smith et al. 2024), but a detailed understanding of how the en-ergetics and observables of outflows depend on resolution remains lacking.
Furthermore, outflows are highly dynamic in nature, propagating at high velocities (⪆ 100 km s −1 ) into highly stratified media where the density and strength of the radiation field steadily decrease away from the galaxy.These quickly-evolving physical conditions and their low-densities make the cooling, chemical and ionic composition of outflows particularly prone to non-equilibrium effects during their propagation, as ionic recombination timescales can significantly differ from the gas cooling time in low-density environments (e.g.Oppenheimer & Schaye 2013, see also Sarkar et al. 2022 for further radiative transfer effects in outflows).
To gain a better understanding of the importance of these physical processes, we create a new suite of galaxy-formation simulations aimed to provide a more robust modelling of galactic outflows.We perform isolated radiation-hydrodynamical simulations of a dwarf galaxy with the AMR code ramses-rtz (Katz 2022).This allows us to (i) obtain an accurate account of non-equilibrium effects by solving the non-equilibrium chemistry of ≥60 ionic species coupled on-the-fly with the spatially-and time-varying radiation field, while (ii) following the development of galaxy-scale outflows selfconsistently powered by stellar feedback within the disc.We further complement this setup by implementing a new AMR refinement strategy in ramses-rtz that explicitly aims to resolve the local gas cooling length (see also Simons et al. 2020 for a similar approach in the CGM), enabling us to better resolve instabilities and mass transfer between gas phases during the outflow propagation (Figure 1, right column).
In this paper, we focus on understanding how gradually improving the resolution in galactic winds impact their energetics, and will dedicate companion papers quantifying how this improved treatment affects their emission lines properties (see also Cameron et al. 2023 for a similar approach with ramses-rtz).We describe the simulation suite and our numerical methods in Section 2. We then show in Section 3.1 that improving outflow resolution enhances the prevalence of the cold and hot gas phases, both of which exhibit systematically increasing energetics (Section 3.2).At fixed feedback budget, this results in a greater than five-fold increase in mass, energy and metal loading factors away from the galaxy (Section 3.3) and a doubling of high-ionization (e.g.O vi) covering fractions (Section 4).We discuss our results and new approach in Section 5 and conclude in Section 6.

NUMERICAL SETUP
We present a suite of radiative hydrodynamical simulations of isolated galaxies, all performed with the same physical modelling but gradually increasing off-the-plane spatial resolution to resolve thermal instabilities in the diffuse outflowing gas.We summarize our numerical galaxy formation model in Section 2.1, in particular how we account for non-equilibrium, ion-by-ion cooling coupled to onthe-fly radiative transfer (see Katz 2022 for a more extensive description of the same setup).Section 2.2 presents the implementation of the new refinement strategy targeting the gas cooling length and Section 2.3 the specific galaxy and simulations of this work.

Chemistry, cooling and galaxy formation physics
We perform radiative hydrodynamical numerical simulations using the adaptive-mesh refinement code ramses-rtz (Katz 2022), built upon the ramses and ramses-rt software (Teyssier 2002;Rosdahl & Teyssier 2015).We solve the gas dynamics using a HLLC Riemann Figure 2), cooling processes are marginally or significantly under-resolved (bottom left, blue).Forcing refinement on  cool (right, here down to 18 pc) resolves the thermodynamics of the outflowing gas much better, leaving only thin interfaces of under-resolved cooling gas outside the disc plane (bottom right).This improved treatment significantly enhances outflow energetics (Section 3) and modifies their ionic structure (Section 4).Accurately capturing the launching and cooling of outflows within the dense ISM requires sub-pc resolution (top, dark red) posing a distinct, complementary challenge for galaxy formation simulations.
solver (Toro et al. 1994), assuming an ideal gas equation of state with adiabatic index of 5/3, while collisionless dynamics of stars and dark matter particles are computed by means of an adaptive particle-mesh method (Guillet & Teyssier 2011).
A key novelty in our simulations is the treatment of ion nonequilibrium thermochemistry, coupled on-the-fly to the radiative transfer.We solve the dynamics of the radiation field at every timestep using the M1 method (Rosdahl et al. 2013;Rosdahl & Teyssier 2015), discretizing its spectrum in eight energy bins from the infrared to the UV (Kimm et al. 2017) and propagating radiation using the reducedspeed-of-light approximation ( reduced = /100).Radiation from local stellar sources is advected using the M1 method, while a uniform and fixed UV background also permeates the simulation box.Both sources contribute to photo-ionization and photo-heating, with the UV background following the tabulation of Haardt & Madau (2012).Gas above  H ≥ 10 −2 cm −3 self-shields exponentially from this UV background (Aubert & Teyssier 2010;Rosdahl & Blaizot 2012).
We then track the non-equilibrium evolution of 11 atomic species Figure 2. Edge-on (top) and face-on (bottom) AMR spatial resolution across a thin slice in the image plane, for the four galaxies with increasingly resolved  cool (left to right) at a time of comparable outflow rates ( = 540 Myr, ≈ 0.05 M ⊙ yr −1 ).The traditional refinement strategy (left) has a strongly layered vertical structure (top left) which under-resolves the off-the-plane cooling length (Figure 1), transitioning quickly from the ISM resolution (18 pc, deep blue) to ≥ 100 pc (white).By contrast, resolution with our new scheme targeting the cooling length naturally tracks the interfaces of expanding superbubbles off the disc plane, following their irregular spatial structure and preventing them from numerically over-mixing.The structure of the dense ISM (bottom, deep blue) only gets visually modified when the additional refinement in the diffuse gas reaches a maximal resolution equal to that in the ISM (right-most).Nonetheless, star-formation conditions between all runs remain consistent with the fiducial setup, even for our most resolved cases (Appendix A). and 64 individual ions (H i − ii, He i − iii, C i − vii, N i − viii, O i − ix, Ne i − vii, Mg i − vii, Si i − vii, S i − vii, Fe i − vii and Ca i) and the formation and evolution of molecular hydrogen (Katz et al. 2017).Number densities of individual ions and species are computed at every simulation timestep by solving a network of coupled equations accounting for recombination, collisional ionization, charge exchange and photo-ionization extracted from the local, time-varying radiation field (Katz 2022).
Having access to individual number densities of ions and the radiation field further allows us to self-consistently compute the out-ofequilibrium cooling rate for the gas.At low temperatures ( ≤ 10 4 K), we explicitly compute the level populations and the fine-structure emission of O i, O iii, C i, C ii, Si i, Si ii, S i, Fe i, Fe ii, Ne ii which dominate the cooling rate at the low metallicities we consider in this work (e.g.Glover & Jappsen 2007).We use tabulated nonequilibrium ion-by-ion cooling tables for our UV background at high temperatures ( ≥ 10 4 K; Oppenheimer & Schaye 2013) and further account for non-equilibrium atomic cooling processes of H and He (Rosdahl et al. 2013) and molecular cooling of H 2 (Katz et al. 2017).
This thermodynamical setup is complemented by an extensive galaxy formation model, including star formation and stellar feedback processes.Stars form in a molecular-cloud-like environment identified using a thermo-turbulent star-formation criterion (Federrath & Klessen 2012;Kimm et al. 2017) and with a Kroupa (2001) initial mass function.Star particles are spawned with minimum initial masses of 1830 M ⊙ .Stellar particles then re-inject energy, mass, momentum, metals and radiation locally, following the mechanical feedback prescription of Kimm et al. (2015) for SNe explosions and the model of Agertz et al. (2013) for stellar winds.Radiation is injected in gas cells surrounding stellar particles assuming a single stellar population of the age and metallicity of the particle and a BPASSv2.2.1 spectral energy distribution (Stanway et al. 2016).The radiation impacts the gas via photo-ionization, radiation pressure, and photo-heating.These simulations model photo-electric and H 2 heating following Katz et al. (2017) and ignore cosmic ray processes.

Refinement strategy on the gas cooling length
All our adaptive-mesh simulations start by using a traditional 'fiducial' quasi-Lagrangian refinement strategy targeting the ISM.To achieve this, we split AMR cells when their dark matter plus baryonic mass exceeds 8,000 M ⊙ for 13 levels in a 150 kpc box, corresponding to a maximum spatial resolution of 18 pc.We also split cells when the local Jeans' length is resolved by less than four cells to obtain a more homogeneous spatial resolution across the ISM.This allows us to better capture ISM turbulence (e.g.Federrath et al. 2011) rather than solely focussing on dense clumps.
Figure 2 illustrates the structure of this fiducial refinement strategy (left-most column), showing an edge-on (top) and face-on (bottom) map of the spatial resolution at a time of comparable outflow rates between the simulations ( = 540 Myr, ≈ 0.05 M ⊙ yr −1 ).With the quasi-Lagrangian approach, resolution is near its allowed maximum within the disc (18 pc, deep blue) but drops rapidly in the vertical direction (≥ 100 pc, white) as the density decreases.
However, the cooling length of gas escaping the galaxy can be much smaller and is either marginally or under-resolved by the quickly degrading resolution (Figure 1).To remedy to this issue, we store the net cooling rate Λ net obtained at each simulation timestep by the non-equilibrium solver and compute for each gas cell.Here,  th and  are the thermal gas pressure and density used to define the isothermal sound speed.The cooling time,  cool , is defined as the internal energy of the gas divided the net cooling rate Λ net , with  the mean molecular weight (mean mass per particle in units of the proton mass), and Λ net in units of erg s −1 cm −3 (see Sutherland & Dopita 1993, eq 52 to 56 for definitions).When computing  cool through Equation 1, we use the self-consistent , Λ net ,  and  computed internally by the non-equilibrium thermochemistry and hydrodynamical solvers.Λ net and  cool can be positive or negative depending on whether the gas is cooling or heating.
We then implement an additional refinement criterion, splitting gas cells in 8 where  cool is positive (i.e.only cooling gas) and unresolved by at least  cool /( cells Δ) < 1, where Δ is the spatial resolution at a given level of the mesh hierarchy.We allow  cells to take different values at different levels along the mesh hierarchy.
The right columns of Figure 2 show the resulting spatial resolutions with this new strategy, resolving  cool with  cells = 8 down to 72, 36 and 18 pc (second, third and fourth columns respectively).The structure of the galactic disc (bottom) and its dense ISM (deep blue) where star-formation occurs is visually similar in all cases, until we reach a cooling length refinement matching the maximal ISM resolution (right-most column).Even in this case however, we show in Appendix A that our new refinement scheme adds resolution in the diffuse gas around the disc, but does not impact the location and densities at which star formation occur and supernovae explode.
Rather, targeting  cool introduces large amount of additional offthe-disc-plane resolution that tracks discrete plumes and shells of outflowing gas as they expand and mix with the inner CGM (top panels of Figure 2).The filamentary structure extending off the galaxy is visibly more extended from left to right, but remains contained by avoiding unnecessary resolution in the hot halo where densities are lower and the cooling lengths larger (white).
Since the ISM is already at the maximum resolution due to the quasi-Lagrangian scheme, the additional refinement criterion on the cooling length allows us to allocate additional resources to the diffuse gas in a controlled and gradual way.This provides us with the ability to cleanly and conveniently decouple the resolution between dense and diffuse gas, and to isolate their relative impact on outflow properties in Section 3.
Furthermore, we also verified that the resolution structure significantly varies with time and naturally follows to the star formation and outflow activity of the galaxy.This provides us with an efficient refinement strategy that focusses on the time-evolving physical quantity of interest rather than a pre-defined, arbitrary volume, and its adaptive nature also enables us to tailor the computational load to the application at hand.One can either model isolated outflows as resolved as the ISM (e.g.our most-resolved case in this study) or use more gradual strategies in the diffuse gas to limit the computational expense (e.g.our intermediate setups).We discuss the computational costs and pros and cons of different foreseen simulation strategies in Section 5.1.

Simulation suite
All simulations in this work (Table 1) follow the evolution of the same, isolated dwarf galaxy, extensively described as G8 in Rosdahl et al. (2015); Kimm et al. (2018); Katz (2022).Briefly, we set up 2 × 10 5 particles to sample a dark matter halo of 10 10 M ⊙ , an initial disc and bulge (masses of 3 × 10 8 M ⊙ and 3 × 10 7 M ⊙ , respectively) leading to a maximum circular velocity of ≈ 30 km s −1 .Disc gas is initialized with an exponential density profile and a weak metallicity gradient peaking at 0.1 Z ⊙ , and is surrounded by a homogeneous (bottom) for our four galaxies with increasingly resolved outflows.Despite the change in numerical setup, all galaxies form a similar amount of stars, ensuring that the total feedback budget is approximately the same across all simulations.Better resolving the cooling length also yields comparable average star formation rates and mass outflow rates close to the disc, providing a controlled study.We discard the first 200 Myr of each simulation to allow the artificial transient phase induced by our initial conditions to dissipate.
( = 10 −6  p cm −3 ) hot, metal-free halo.We initialize all ions in their ground state in the disc and in their most ionized state in the hot halo.
We perform four main simulations of this galaxy, progressively increasing the resolution target of the cooling length to 72, 36, and 18 pc eventually matching the ISM resolution.We follow the evolution of each galaxy for 750 Myr (i.e.∼ 5 disc crossing times) and save snapshots at least every 3 Myr.Figure 3 shows the evolution of the stellar mass formed over the course of each simulation (top), their star-formation rates averaged over 10 Myr (middle) and mass outflow rates close to the disc of the galaxy (bottom, measured in a 0.1 kpc-thick slab placed at || = 1 kpc and of cylindrical radial extent  = 4 kpc; see Figure 4 for a visualization).
Our initial conditions drive an initial starburst which recedes after ≈ 100 Myr, and it takes an additional ≈ 100 Myr for the galactic disc and outflow to settle after a large enhancement in mass outflow rate (grey shaded region in Figure 3).We thus start our analysis at 200 Myr, ensuring at least 130 snapshots of analysis when deriving time-averaged properties (Section 3).We verified that our conclusions are unchanged if starting at 250 Myr or 300 Myr.
All simulations form a similar amount of stellar mass, ending Table 1.Summary of the simulations presented in this work (first column).
We simulate the same initial conditions with a maximum resolution in the ISM of 18 pc, refining the outflows to an incrementally higher spatial resolution (second column).Two of our setups are simulated three times with different random seeds (first and third lines) to assess stochasticity in our results (Appendix B).We also provide a controlled, time-limited test of the computational cost of each setup (third column) discussed in Section 5.1.up within 15 per cent of each other.This is key to compare their outflows, as it guarantees that the overall energy, momentum and metal budget from stellar feedback is nearly unchanged between different resolution runs.It also provides further evidence that the physics of star formation in our simulations is not affected by the additional resolution in the diffuse gas (Appendix A).However, the ability to drive galactic outflows also depends on the coupling between feedback and the surrounding gas, which in turn depends on the clustering of SNe and the density of the medium in which they explode.This is subject to stochasticity due to the probabilistic nature of our modelling of star formation (see Genel et al. 2019;Keller et al. 2019 for further discussion).To quantify the magnitude of this scatter, we perform three simulations of our reference galaxy and three with the cooling length resolved down to 36 pc, solely varying the random seed for star formation in each of them.We show in Appendix A that star formation and feedback conditions within the ISM are comparable between both stochastic resimulations of a given setup and different resolutions of the cooling length.We further show in Appendix B that run-to-run stochasticity leads to different but statistically compatible outflow rates close to the disc, and that the induced scatter remains subdominant compared to differences in outflow properties and rates resulting from an increase in resolution when measured further away from the galaxy (Section 3).This therefore confirms the causal association of differences in outflows with a more accurate treatment of their thermodynamics.

ENERGETICS OF BETTER-RESOLVED OUTFLOWS
We start by visualizing the impact of better resolving the cooling outflow from the galaxy.Figure 4 shows thin slices of gas density (top) and temperature (bottom) in a side-on view at the same time as in Figure 2. As the cooling length of the diffuse gas is increasingly resolved (left to right), the thermodynamical structure of the outflow is visibly affected.Gas goes from a smooth and warm medium (left) to an increasingly structured distribution, with both colder and higherdensity shells and hotter and lower density cavities (right).As we will see in Section 3.1, increased density contrasts and a more multiphase structure are a systematic consequence of the improved numerical treatment of the diffuse gas and extends beyond the specific time of Figure 4. Furthermore, in the bottom panels, we overlay velocity streamlines to visualize changes in gas kinematics.At this time, the streamlines also showcase greater opening angles with increasing resolution in the diffuse gas.Increased opening angles can boost the acceleration of galactic winds (e.g.discussion in Martizzi et al. 2016), although Appendix C shows considerable variability over time in A better resolved cooling length (purple to green) leads to longer tails towards both denser and more diffuse gas, and both colder and hotter gas.This points towards an increasingly multiphase structure of galactic outflows with increasing resolution (Figure 6).
opening angles that makes systemic trends with numerical resolution difficult to establish.We show in Section 3.2 that better-resolved galactic outflows are indeed systematically faster moving and more energetic, significantly boosting outflow loading factors through the slab and shell visualized in Figure 4 as a result (Section 3.3).

Increasingly multiphase outflows
We first quantify the thermodynamical structure of our outflows over time.To this end, we extract a cylinder of gas above and below the disc plane ( ≤ 4.0 kpc and 0.5 ≤ || ≤ 5 kpc) and select outflowing gas with (  ) = ()1 using || > 1.0 kpc, || ≤ 4 kpc, || ≤ 7 kpc or  ≤ 5.0 kpc..At each simulation snapshot, we then compute the mass-weighted two-dimensional histograms in pressure-density and temperature-density (150x150 pixels) of this gas.We stack these histograms pixel-by-pixel over ≥ 500 Myr of time evolution, and normalize by the total gas mass to obtain the one-and two-dimensional PDFs in Figure 5 and Figure 6.Despite the change in resolution, all runs have similar total masses of outflowing gas (varying by less than 20 per cent, top of each panel in Figure 6), ensuring that changes in the PDFs are largely unaffected by changes in overall normalizations.
Starting from the one-dimensional distributions for our fiducial case (Figure 5, purple), we see that the majority of the outflowing gas mass is at low densities (∼ 10 −3 −10 −1  p cm −3 ) and warm temperatures (∼ 10 4 K), as expected from diffuse gas in photo-ionization equilibrium with the surrounding UV background.A subdominant fraction of the gas is dense enough to self-shield efficiently against the UV background ( ≥ 10 −1  p cm −3 in our model) and cool below 10 4 K. SN-heated gas is also visible as a long tail towards high temperatures (≥ 10 5 K), while the enhancement of gas with  ≈ 8 × 10 4 K maps to a local minimum of the cooling function at the low metallicities considered in this work (≤ 0.3 Z ⊙ ; e.g.Bialy & Sternberg 2019;Katz 2022;Kim et al. 2023).
Improving the numerical treatment of thermal instabilities (purple to green) does not modify this broad picture but introduces systematic trends for the tails of these distributions.The high-density and low-temperature tails incrementally reach higher densities and lower temperatures, confirming the expectation that condensation and cooling in the diffuse gas are numerically suppressed in our fiducial setup (Figure 1).However, the improving resolution also generates a longer tail towards higher temperatures and lower densities, highlighting that the hot, diffuse phase is becoming increasingly significant.Since the total gas mass is roughly conserved between simulations, this is mainly modifying the relative fractions between gas phases, and thus the multiphase structure of our galactic outflows.
To better visualize these changes, Figure 6 shows the temperaturedensity (top) and pressure-density (bottom) two-dimensional phase diagrams of the same gas as in Figure 5, recovering the broad features identified in the one-dimensional PDFs (top, annotated).
With increasing resolution (towards the right), the low-density gas ( ≤ 10 −2  p cm −3 ) is better sampled, generating a larger and smoother scatter in pressure and temperatures for a given density.This low-density gas also has a visibly cooler tail (top right), which we check arises from better-resolved ionized cavities following SNe explosions.These cavities cool from adiabatic expansion in a regime of inefficient recombination, preventing rapid heating from the UV background (see also McQuinn 2016 for a similar argument driving the temperature-density relation of the IGM).The dense gas ( ≥ 10 −1  p cm −3 ) cools more efficiently below 10 4 K with increasing resolution (top right), generating a more extended scatter in temperatures and starting to generate two distinct pressure tracks.
Our results corroborate the physical picture in which the cooling length of diffuse outflowing gas from galaxies is under-resolved with traditional simulations (Figure 1 and Section 1), hampering the development of multiphase galactic outflows.Improving resolution in the outflowing gas significantly helps alleviate this issue, yielding a more prominent colder and denser phase, but also a more prominent hotter and diffuse phase as we avoid over-mixing and over-cooling .In all cases, most gas mass is warm and in photo-ionization equilibrium with the UV background (≈ 10 4 K).However, an increasing resolution yields a better-sampled diffuse phase with an increased scatter in temperatures and pressures at a given density, and allows denser and colder gas at high densities.This increasingly multiphase structure in turn affects how the thermo-kinematics of outflows as they propagate from close to the disc (Figure 7) into the inner CGM (Figure 8).
SN-driven expanding superbubbles.Both of these components are particularly important for galactic-scale feedback, as the cold-towarm phase is expected to be significantly mass loaded, while the hot phase deposit energy in the surrounding diffuse CGM.We quantify the change in energetics of each of these gas phases in the next Section to gain better insights into these aspects.Before turning to this, we briefly highlight that the one-or twodimensional PDFs do not show signs of numerical convergence with increasing resolution and that the ordering of simulations with resolution is not linear (e.g.blue versus green in Figure 5).This hints that, despite improving the situation, we are yet to fully resolve the multiphase structure of the outflowing gas and that evolution at a given resolution possesses stochasticity due to the specific star formation history at hand.We quantify this stochasticity in Appendix B, showing that it is a sub-dominant effect compared to the broad trends we identify when increasing resolution, and discuss remaining numerical uncertainties and convergence in Section 5.2.

Kinematics and energetics of better resolved outflows
To gain physical insights into these differences in outflow energetics, we show in Figure 7 the two-dimensional distributions of temperature and vertical velocities of outflowing gas close to the disc (a 0.1-kpc slab at || = 1kpc with (  ) = ()).The histograms are weighted by gas mass and total gas specific energy (defined as the sum of the kinetic and internal specific energies,  gas = 1/2  2 +  = 1/2  2 +  th //( − 1)) in the top and bottom panels, respectively.Again, each panel shows the normalized stack over the time evolution of each galaxy, with increasingly resolved cooling length from left to right.
From Figure 7, we recover that most of the mass-weighted out-flow gas (top panels) is warm (≈ 10 4 K), travelling at relatively low velocities (  ≤ 10 km s −1 ) and is unlikely to escape the gravitational potential of the galaxy ( circ as grey dotted line).This picture remains qualitatively similar as the cooling length is increasingly resolved (left to right), although the improved resolution leads to a tail of ≤ 10 4  gas (top, left panel).This confirms that, close to the disc, the ejection of mass by SNe-driven outflows is not significantly modified by adding additional resolution in the diffuse gas (recall Figure 3 and see also Appendix A).
In stark contrast, the improved numerical treatment of the hot phase at higher resolution yields significant differences when weighting the same histograms by gas specific energy (bottom panels).Focussing first on the fiducial case (left), the hot, fast-moving and energycarrying gas phase of the outflow materializes as a track extending towards high temperatures and velocities, broadly following an iso-Mach slope (dashed lines, where M = /  with   = √︁  th / the isothermal sound speed).This 'plume' arises from the superposition of the many SN-driven superbubbles expanding out of the disc plane (see also Kim et al. 2020;Andersson et al. 2023;Steinwandel et al. 2024).This gas is strongly supersonic, and dominated by its kinetic energy (solid line shows the equality between gas kinetic and internal energy).Gas sampling the radiative regions behind shocks (hot and subsonic, towards the top left) is present, but poorly sampled in the fiducial case.
As we increasingly resolve the cooling length in diffuse gas (left to right), the plume of hot and fast-moving gas extends towards faster velocities and higher temperatures (respective upper-right corners of each panel).The hot phase is thus launched more energetically close to the disc.Furthermore, a significant hot, but subsonic, phase becomes ever more apparent with each level of increased resolution in the cooling length.This is to be expected as, by design, our approach concentrate resolution in rapidly cooling regions such as the Each level of increased resolution also makes the cloud of hot and subsonic gas more prevalent, indicating that we increasingly resolve the sonic transition of the hot phase.This is key to accurately capturing the acceleration of the hot gas and its subsequent energetics as it escapes into the CGM (Figure 8).  4 for a visualization).At these outer radii, outflow mass is hotter and faster moving ( ≥ 10 km s −1 ) than at |  | = 1 kpc.The improved numerical treatment of the thermo-kinematics of the gas significantly boosts its energetics, with outflow mass systematically travelling at faster velocities (top panels) and a more energy-loaded hot-fast phase (bottom panels).The combination of these factors translate into large boosts in outflow loading factors (Figure 9).post-shock, radiative layers that generate hot and fairly slow-moving (  ≈ 50 km s −1 ) gas.Combining these aspects, we conclude that our improved resolution scheme drastically improves the treatment of the thermo-kinematics of the hot phase by better capturing the expansion and cooling around shocked regions and keeping gas hotter for longer by avoiding spurious mixing and numerically-driven rapid cooling.

Temperature (K)
We can expect such improvements at  = 1 kpc to be reflected in the energetics of the wind at larger radii, as better capturing the transition of the hot phase from sonic to supersonic is key to capture the global acceleration of galactic outflows (e.g.Chevalier & Clegg 1985;Smith et al. 2024).To visualize this, we repeat the same computations as for Figure 7, but now selecting gas at  = 5 kpc (≈ 0.2  200 ) and plotting   - to quantify differences in outflow thermo-kinematics after propagation into the inner CGM.In Figure 8, we only include gas within a biconical shell with a broad opening angle  = 45 • to capture the opening of the velocity streamlines.We show in Appendix D that this selection removes gas from polar angles close to the plane of the disc that is often dominated by the hot hydrostatic halo remaining from the initial conditions of our galaxy (rather than the outflowing gas of interest).
At  = 5 kpc, the mass-weighted outflow (Figure 8, top panels) is hotter and faster moving than at  = 1 kpc, as expected if the gas is accelerating as it escapes the galaxy.Furthermore, marginalizing the mass-weighted diagrams of Figure 7 and 8 along the y-axis, we then compute the mass-weighted velocity PDF.We find that both these PDFs shifts towards larger velocities with increasing resolution, but the increase is relatively mild at  = 1 kpc (medians going from 6.7 to 8.5 to 8.5 to 9.2 km s −1 with each added resolution level, respectively) and more significant at 5 kpc (going from 39.9 to 63.4 to 46.5 to 50.3 km s −1 ).In all cases, the outflowing gas significantly accelerates between  = 1 kpc and  = 5 kpc, but this acceleration is amplified with higher resolution outflows.
Weighting the outflow by energy confirms this picture (bottom panels in Figure 8).Even after propagating to  = 5 kpc, the hot phase remains better resolved with our improved scheme, extending towards hotter temperatures and faster velocities as expected from better capturing hot-to-warm interfaces (see e.g.Kim & Ostriker 2018 for a discussion).At our highest resolution (right column), we even recover a small hot, subsonic component (although much less prominent than at  = 1 kpc) that is entirely absent with the fiducial resolution scheme.This highlights that we now capture the thermodynamics of shocks and the sonic transitions of specific superbubbles in the wind out to large radii.As we will see in the next Section, this results in significant boosts to the outflow loading factors with resolution.

Consequences on outflow loading factors
We now turn to quantifying the impact of the improved numerical treatment of the outflowing gas on its loading factors.We compute, at each saved simulated time, the mass, energy and metal loading factors as Here, S is the interface with thickness ΔS through which the outflow rate is defined (respectively, a 0.1-kpc thick horizontal slab placed 1 kpc above and below the disc, and a 2-kpc thick shell placed at  = 5 kpc and polar angle  ≤ 45 • ; see Figure 4).  ,   ,   and   are the mass, specific internal energies, velocities and total metal mass fraction of individual gas cells within S. We use   =  , in Equation 2 for the slab and   =   , for the shell to account for the opening of the streamlines (Figure 4).We normalize the outflow loading factors with the star formation rate of the galaxy averaged over the last 10 Myr2 , the average specific energy injected by SNe for a Kroupa (2002) mass function,  SN = 5.2 × 10 5 km 2 s −2 (Kim & Ostriker 2018) and the average mass-weighted metal mass fraction within the disc  disc (geometrically defined as gas with  ≤ 4.0 kpc and || ≤ 0.5 kpc).Relative trends between outflow loading factors are unchanged if measuring outflow rates at || = 0.5 kpc or using other radial definitions (see Appendix D).
Figure 9 shows the mass, energy and metal outflow rates (top, middle and bottom panels respectively) at 1 and 5 kpc (left and right panels respectively).Violin plots showcase distributions across the time evolution of each simulation, highlighting their medians (diamonds) and 16-84 confidence intervals (black lines).
Outflows close to the disc (left panels) are significantly mass and metal loaded (  ,   ≥ 1) but not energy loaded (  ≪ 1).Enhancing off-the-disc-plane resolution leads to a mild but systematic evolution, with a 50 per cent enhancement in the median mass loading factor between the fiducial and the most resolved case, a factor of two growth in the energy loading factor, and a 20 per cent increase in the median metal loading factor.These limited increases in mass and metal loading are consistent with the relatively unchanged mass-weighted PDF in Figure 7 (top panels), while the increase in energy loading factor is consistent with the hotter hot phase and better captured transition of the sonic-to-subsonic transition close to the disc (Figure 7, bottom panels).However, except for the energy loading factor, the significance of these increases remains difficult to interpret, as stochasticity due to varying star formation activity at different times (black lines show the 16-84 confidence intervals) as well as the overall noise in the evolution of a numerical setup at a given resolution (Appendix B) can generate shifts of a similar magnitude.
In contrast, differences further away from the galaxy (right panels) are far more significant, with a two-fold increase in mass loading (median   = 0.20, 0.34, 0.43, 0.47 from purple to green) and metal loading factors (median   = 0.46, 0.72, 0.84, 0.98), and a five-fold increase in energy loading factor (median   = 0.0034, 0.0085, 0.0090, 0.017).Importantly, these shifts with enhanced resolution are close to exceeding the 16-84 confidence interval span of the fiducial case and are robust against run-to-run stochasticity (Appendix B).
Our study shows that better resolving the diffuse outflowing gas provides a significant boost to their energetics and, in turn, to their ability to regulate galactic star formation.When escaping into the CGM, our increasingly-resolved outflows go from being marginally ( ≤ 1) to significantly ( ≥ 1) mass and metal loaded.This ability to maintain the high mass and metal loadings acquired in the ISM over super-galactic scales could have strong consequences for galactic outflows' ability to efficiently transport mass outwards and pollute the CGM and IGM with metals.Similarly, the growth towards well Better resolving the dynamics of multiphase outflows from the low resolution fiducial case (purple) to the highest resolution case (green) systematically increases all loading factors.Close to the disc (left), increases are limited for the mass and metal loading factors, but grow to a two-fold increase at larger radii (right).The better-captured propagation of the hot phase significantly enhances the energy loading factors through both interfaces, culminating in a 5× boost at 5 kpc.
energy-loaded outflows at 5 kpc would naturally help pressurize the inner CGM and prevent further gas accretion onto the galaxy.
The full reach of these results remains to be determined in a cosmological context, however.Our outflows currently propagate in a low-density, static hot halo surrounding our idealized dwarf galaxy, rather than the layered, dynamic and inflowing CGM of a cosmological galaxy.However, their significance strongly motivates future studies quantifying how the gas cycle of galaxies is affected when outflows and inflows are better resolved.

IONIC STRUCTURE OF BETTER-RESOLVED OUTFLOWS
A key novel aspect of our simulations is to track the non-equilibrium abundances of ions, providing us with the opportunity to readily study the ionic structure of our outflows and eventually link to emission and absorption observables.We now provide preliminary insights on the observational consequences of our findings, and will study specific observational diagnoses to recover gas properties from outflow observations in dedicated companion papers (see Section 6 for further discussion).
Figure 10 shows the time-averaged covering fractions of five ions with increasing ionization potentials.To derive those, we generate side-on images of the galaxy at each time output (8 kpc wide and deep, spatial resolution of 20 pc) masking the galactic disc (|| < 0.5 kpc) and compute histograms of line-of-sight column densities for each ion.We then aggregate the histograms over time, renormalize them to obtain the stacked probability distribution functions, and take the cumulative sums shown in Figure 10.

OVII
Figure 10.Time-averaged covering fractions of selected ions with increasing ionization potentials (left to right panels respectively) when looking at simulated galaxies edge-on in runs with increasing off-the-plane resolution (from purple for the fiducial case to green for the most resolved: see text for detail).A better resolved warm-to-cold phase leads to an increase in Hiand Mg ii, particularly at low column densities (≤ 10 17 cm −2 and ≤ 10 12 cm −2 ).A better resolved hot phase strongly enhances the covering fractions of high-ionization ions, notably in C ivand O vi, two key ions to probe CGM gas in absorption.
tween our two extreme resolutions by 5 and 8 per cent, respectively.This is primarily driven by an increase at high column densities (e.g. ( Hi ≥ 10 18 cm −2 ) increases by 15 per cent), which is best explained by the growth of a colder and denser gas phase as the out-of-the-plane cooling length is better resolved (Figure 6).The relatively limited evolution in the covering fractions of these lowionization ions with resolution is encouraging, although the lack of convergence highlights that we are not yet resolving the characteristic scales regulating the dense, cold phase (see e.g.van de Voort et al. 2019; Gronke & Oh 2020 and Section 5.2 for further discussion).
For ions tracking warmer, more diffuse and ionized gas (e.g.C iv, middle panel), we find a reduction in high-column density sightlines and a growth at low column densities with increasing resolution.We attribute the latter to the larger scatter of temperatures in the diffuse phase (Figure 6) yielding a more extended scatter in diffuse C iv sightlines.In contrast, higher-densities sightlines hosting C iv at low resolution have likely recombined and cooled to lower-ionization ions when the cooling length is better resolved, explaining the reduction.
A much larger effect is visible when inspecting the covering fractions of higher ionization ions (right, O vi, O vii), with  ( O vi ≥ 10 13 cm −2 ) going from to 0.148 to 0.160 to 0.169 to 0.246 with increasing resolution (0.306 to 0.321 to 0.327 to 0.437 for O vii).These 66 and 43 per cent increases materialize the increasing energetics of our better-resolved outflows, with the systematically hotter gas promoting higher-ionization ions.
The trend with resolution is however non-linear and difficult to interpret quantitatively.Lower-resolution runs (violet and blue) show mild evolution until a sudden increase in the covering fraction occurs with our most-resolved run.We checked that this is in part due to stochasticity in the given realization of the star formation and outflow history.For example, our two other realizations of the outflow at 36 pc (Appendix B) show 74 and 24 per cent increase in  ( O vi ≥ 10 13 cm −2 ) compared to our fiducial simulation, while our other two random realizations of the fiducial setup vary by as much 27 per cent.Despite these quantitative uncertainties, we find that runs with improved treatment of thermal instabilities in the diffuse gas always have elevated covering fractions of O vi and O vii.
Our results show that improving the treatment of thermal instabilities in the diffuse gas provides a significant avenue to increasing the covering fraction of high ionization ions.Probing the thermodynamics of the gas surrounding dwarf galaxies using ionic absorption lines is already possible (e.g.Bordoloi et al. 2014;Johnson et al. 2017;Zheng et al. 2019Zheng et al. , 2020Zheng et al. , 2024;;Qu & Bregman 2022), although we stress that our results should be compared with great care (if at all) to such data.Our idealized simulations lack a self-consistent large-scale CGM surrounding the dwarf, strongly limiting comparisons at large impact parameters and against statistical samples of dwarf galaxies.We verified that our results are qualitatively unchanged if including slightly higher impact parameters (using 10 kpc-wide images rather than 8), but refrain from detailed comparisons with data.
Nonetheless, our results stress the importance of improving out-ofthe-plane resolution to properly capture the ionic structure of outflows and constrain galaxy formation models from interpreting their observables in emission and absorption.These results further motivate applications of our setup and refinement strategy in a cosmological context to provide more robust comparisons against observational measurements at larger impact parameters, well into the CGM.

The numerical cost of refining on the cooling length
To establish that better resolving the multiphase structure of galactic outflows boosts their energetics, we rely on adding resolution in the diffuse thermally-unstable gas.This carries an associated numerical cost that potentially limits the practical scope of our approach and the applicability of our refinement strategy.We discuss in this section the cost of refining on the cooling length and strategies we envision to mitigate it.
Simulations in our main suite were performed with a different number of cores across different machines.Their total costs thus mix the increased computational load with different intrinsic hardware and parallel efficiencies, making a direct comparison difficult.To provide a cleaner test, we start from our most refined setup at the time shown in Figure 1 and 2 ( = 540 Myr; i.e. neither particularly quiescent nor outflowing) and restart four runs for 6 wall-time hours on 768 cores of the same machine, gradually degrading the cooling length refinement target to 36 pc, 72 pc, and none.Table 1 (right column) reports the median timestep cost (in CPU wall-clock time per simulated Myr) for each of these tests, which we stress can not be extrapolated to estimate the total simulation cost as the timestep strongly varies depending on the outflowing conditions.
As expected, the median cost of a timestep is rising as we increasingly resolve the cooling length, driven by the accordingly increasing number of grid cells (Figure 2).However, this cost only increases by 2 and 28 per cent for our intermediate setups that resolve the cooling length to 72 and 36 pc respectively.It then sharply grows to a 3-fold increase in the most refined run that matches the cooling length with the ISM resolution at 18 pc.This non-linear increase reflects the fact that the computational effort is dominated by the densest, most refined cells that have more restrictive Courant-Friedrichs-Lewy conditions and higher convergence costs for the chemistry and radiative transfer solvers.In other words, intermediate setups resolving the diffuse gas twice or four times worse than the ISM are not free, but their computational costs can be kept under control as adding many coarser cells remains a sub-dominant cost compared to computations in the fine ISM gas cells.Given the already large gains in outflow loading factors we observe at 72 and 36 pc (Figure 9), we foresee these intermediate setups as the most attractive and tractable options for cosmological setups.
In fact, cosmological simulations of galaxies with such an AMR strategy refining on the cooling length have already been achieved (Simons et al. 2020;Lochhaas et al. 2021Lochhaas et al. , 2023)), although at a comparatively much lower spatial resolution in both the ISM and the cooling length (≈ 200 pc) than here.This lower resolution arises from targeting larger galaxies than the dwarf considered here, and from fundamental differences between the isolated setup of this work and a live cosmological environment.The presence of dense filaments and inflows, of merging galaxies with their stripped tails of gas, and of an already multiphase and turbulent CGM surrounding the central galaxy all add gas highly susceptible to thermal instabilities.This would likely trigger large amounts of additional refinement and sharply increase the computational and memory costs of the simulation.
Although field-testing is required to pinpoint the achievable cooling-length resolution in a given cosmological setup with ramsesrtz, we remain confident of its overall tractability due to the overwhelmingly dominant cost of the ISM gas.Furthermore, reducing the computational expense is possible by relaxing our requirement that the cooling length is resolved by at least  cells = 8, although further exploration is needed to understand the convergence properties of our simulations with  cells .Another mitigation strategy could also be to release additional resolution levels gradually over cosmic time for the cooling length as is often done for the ISM (e.g.Snaith et al. 2018).Or more simply to restart a simulation for a limited time period around a specific event such as a starburst to focus all the effort on one maximally-resolved galactic outflow in a self-consistent environment.

Numerical uncertainties in resolving the multiphase structure of galactic outflows
Our results show that incrementally increasing resolution in the diffuse, outflowing gas boosts both the mass and energy loading factors of outflows as they escape the galaxy.In Section 3, we tie this to a better-resolved multiphase structure.In this Section, we discuss the physical and numerical mechanisms responsible for this effect, and the remaining uncertainties associated with modelling the hydrodynamics of galactic outflows.Multiple physical processes are responsible for setting the multiphase structure of gas on small scales, ranging from shocks, hydrodynamical instabilities and cooling and heating processes (Section 1 and Faucher-Giguere & Oh 2023 for a review).The primary target of this study is to better capture the poorly-resolved cooling length in the diffuse gas (Figure 1).This can lead to an artificial suppression of multiphase structures (e.g.discussion in Hummels et al. 2019).We tackle the issue using the refinement strategy described in Sec-tion 2.2, but a direct consequence of increasing resolution in the diffuse gas is to also help mitigate other numerical effects.
In particular, numerical diffusion of hydrodynamical quantities due to the advection of a fluid in motion with respect to a fixed grid (e.g.Robertson et al. 2010;Teyssier 2015) is particularly important in the context of fast outflowing gas (≥ 100 km s −1 ).Advection errors are demonstrably reduced with increasing resolution, and our most resolved simulations would thus suffer less from advection-driven numerical mixing between gas phases.Similarly, the heating, lifetime and propagation of shocks are also better captured with increasing resolution (e.g.Skillman et al. 2008;Creasey et al. 2011;Vazza et al. 2011) and typically under-resolved in the diffuse gas (e.g.Bennett & Sĳacki 2020).A better treatment of shocks then comes as a byproduct of refining on the cooling length (e.g. Figure 4), improving their ability to thermalize gas at higher temperatures in the lowest density regions which are better sampled.Furthermore, resolving the transition between subsonic and supersonic gas in the outflow is key to accurately capture the acceleration of the hot phase (e.g.Chevalier & Clegg 1985;Fielding & Bryan 2022).This can be under-resolved with traditional schemes (Smith et al. 2024), and is improved in our more resolved runs (Figure 7).Lastly, hydrodynamical turbulence is heavily suppressed by insufficient resolution, and sets the characteristic growth timescales of the cold phase (Gronke & Oh 2018;Tan et al. 2021;Gronke et al. 2022).The traditional quasi-Lagrangian strategy erases small-scale information in the diffuse gas and likely underestimates the amount of gas turbulence (see also Bennett & Sĳacki 2020).The added resolution visibly improves the level of turbulence in the outflowing gas (e.g. in this movie), despite our refinement strategy still suppressing turbulence due to refinement/derefinement noise (Teyssier 2015, see also Martin-Alvarez et al. 2022 for other strategies to mitigate this issue in a galaxy formation context).
Pinpointing the respective importance of each of these coupled numerical effects remains difficult in our current setup.Disentangling them could be possible by implementing complementary refinement strategies that target one effect at a time (e.g.shock refinement; Bennett & Sĳacki 2020) and are used in turn.But we stress that the numerical gains associated with targeting the cooling length of the diffuse gas are multifold, and extend well beyond better-captured thermal instabilities.
Despite these gains however, our simulations do not exhibit signs of numerical convergence in gas temperature, density or velocity PDFs (Figures 6 and 8), highlighting that we are still under-resolving relevant physical processes.This should be expected as resolving their characteristic length scales of each individual microscopic process is required to obtain numerically-converged results and can require sub-pc resolution in the diffuse gas (e.g.Koyama & Inutsuka 2004;Kim & Kim 2013).However, the interaction of many processes, each with their own preferred regime and characteristic scales, somewhat blurs this picture, and quantities integrated across the plasma might converge without ever resolving microscopic scales (Tan et al. 2021).
In fact, despite the clear growth of outflow loading factors at 5 kpc with increasing resolution (Figure 9), this growth slows as the resolution keeps improving.This hints that convergence in such integrated quantities could be achieved, in particular in a statistical sense, given the stochasticity associated with star formation activity.More quantitative statements, however, require a wider exploration of numerical resolutions and models.
To summarize, the large boost in outflow loading factors induced by better-resolving the structure of the wind (Figure 9) strongly motivates extending our analysis.In particular, future studies quantifying the respective roles of the physical and numerical effects discussed above will be key to gain an understanding on the convergence and robustness of our modelling of SN-driven winds, and of their predicted efficiency at regulating galactic star formation.We stress that such studies would need to be complemented by parallel, complementary explorations of the uncertain launching physics of the wind within the ISM.This aspect remains unexplored in this work, but different star formation models, feedback channels and stellar evolution processes can also generate order-of-magnitude variations in dwarf galaxy outflow loading factors (e.g.Agertz et al. 2020;Emerick et al. 2020;Smith et al. 2021;Andersson et al. 2023;Steinwandel et al. 2023 for recent studies; Naab & Ostriker 2017 for a review).

CONCLUSION
We perform radiative-hydrodynamical numerical simulations of an isolated dwarf galaxy, improving the thermodynamical modelling of the outflowing gas escaping the galaxy.Each simulation drives selfconsistent galactic outflows from SNe and photo-ionization feedback within the ISM resolved at 18 pc, and accounts for the nonequilibrium chemistry and cooling of gas coupled to the radiative transfer using the AMR ramses-rtz code (Katz 2022).
To remedy the under-resolved cooling length in the diffuse outflow when using a traditional quasi-Lagrangian refinement strategy (Figure 1), we implement a new refinement strategy in ramses-rtz targeting the local gas cooling length computed on-the-fly from the non-equilibrium, line-by-line cooling rate.We then perform additional simulations of our dwarf galaxy incrementally resolving the cooling length down to 72, 36 and 18 pc.
The additional refinement criterion on the cooling length improves the resolution in the diffuse out-of-the-galaxy-plane gas (Figure 2), without significantly modifying the average star formation rate, the amount of stellar mass formed (Figure 3), or the star formation and supernova feedback conditions in the dense, ISM gas (Appendix A).This provides us with a controlled study, where the launching mechanics of the outflows within the galaxy are comparable but their propagation as they escape the galaxy is better resolved.
With increasing resolution, outflows become increasingly multiphase, exhibiting both a systematically higher fraction of colder and denser gas, and a more diffuse phase featuring an increased scatter in temperature (Figure 5 and 6).Furthermore, our improved numerical scheme better captures the propagation of the hot and energetic phase of the outflow.Preventing artificial mixing between gas phases leads to systematically hotter gas close to the disc, as shock heating become ever more efficient.Furthermore, better sampling the diffuse gas also allows us to resolve the subsonic, hot regions behind shocks and the key transition between supersonic and subsonic gas (Figure 7).These improvements close to the disc then extend to increased energetics at larger radii, with the hot phase staying systematically hotter at large radii (Figure 8).
Over ≈ 500 Myr of evolution, the combination of these factors leads to a five-fold increase in the average energy loading factor of the galaxy, and two-fold increases in the average mass and metal loading factors in the inner CGM ( = 5 kpc, Figure 9).These boosts occur without modifications to the internal feedback budget or modelling and are robust to stochastic realizations of the star formation history (Appendix B).Rather, they stem from better resolving the hydrodynamical and thermodynamical processes as the outflow propagates into the CGM (Section 5.2).The significance of this boost in outflow energetics without new feedback mechanisms strongly motivates future studies extending the analysis to a cosmological context.Quantifying how such outflows interact with a self-consistent, already-multiphase cosmological CGM, rather than the isolated hot halo of this study, will be key to understanding their improved efficiency at regulating galactic star formation over cosmic time.
Furthermore, a key novelty of our simulations is to track the nonequilibrium abundances of ≥ 60 ions, allowing us to robustly capture the ionic structure of the diffuse gas (Figure 10).Covering fraction of low-ionization ions (e.g.Hi) exhibit a limited increase, primarily driven by a more prominent dense and colder phase better (e.g. Hi ≥ 10 18 cm −2 increasing by 15 per cent).Higher-ionization ions (O vi, O vii) show larger enhancements following the hotter temperatures of our outflows with increased resolution, with both  O vi ≥ 10 13 cm −2 and  O vii ≥ 10 13 cm −2 increasing by ≈ 50 per cent between our extreme resolutions.
Our results demonstrate the importance of additional resolution in the diffuse gas to capture the ionic structure of outflows, and interpret their observable signatures in emission and absorption.Dedicated companion papers building upon the improved treatment of outflows presented here will quantify the robustness of inferences of outflow properties from observational diagnoses (see also this attached movie).In particular, we will post-process our simulations with resonant-line radiative transfer using the rascas code (Michel-Dansac et al. 2020) to obtain mock Mg ii, Si ii and Fe ii spectra of our outflows.This will allow us to assess the accuracy with which such observations (e.g.Rupke et al. 2005;Martin & Bouché 2009;Rubin et al. 2011;Martin et al. 2013) can recover outflow properties and loading factors (H.Katz et al. in prep).We also plan to use our simulations to assess potential biases in recovering spatially-resolved outflow metallicities from emission line maps (Cameron et al. 2021, A. Cameron et al. in prep).
Finally, thermally-unstable and multiphase gas are common features across galaxy formation and astrophysics.The approach presented in this study provides a highly modular AMR strategy to better resolve the cooling length of the diffuse gas, with a computational cost that can be adapted to galaxy formation applications while taking into account the resources at hand (see Section 5.1 for a discussion).We foresee a wide range of potential applications where the combination of this strategy with the non-equilibrium, ion-by-ion cooling of ramses-rtz could provide significant modelling improvements.Namely, to help resolve the multiphase structure of the diffuse, low-density CGM surrounding galaxies (e.g.Hummels et al. 2019;Peeples et al. 2019;Suresh et al. 2019;van de Voort et al. 2019;Lochhaas et al. 2021Lochhaas et al. , 2023)), to better capture the potential shattering of cosmological filaments and sheets (Mandelker et al. 2020(Mandelker et al. , 2021)), or to provide an efficient way to improve the numerical treatment of cosmological ram-pressure striping (e.g.Simons et al. 2020) and better understand the cooling, star-forming tails of jellyfish galaxies (e.g.Tonnesen & Bryan 2012;Tonnesen 2019;Lee et al. 2022).

APPENDIX A: STAR FORMATION STATISTICS
We check in this Appendix that the star formation properties of our simulations are unaffected by the additional refinement scheme on the cooling length, ensuring that the stellar distribution and its associated feedback are consistent between runs.
Figure 3 shows that the total amount of stars formed over the course of the simulation is similar to within 15 per cent between runs at different resolutions.We also checked that the stochastic re-simulations at the same resolution presented in Appendix B verify this condition as well.This guarantees that the total, integrated feedback budget over the evolution of our galaxies is largely unchanged, but not that its instantaneous coupling and efficiency is unaffected.
A formal quantification of this would require storing the properties of each individual SN to obtain their three-dimensional clustering and explosion conditions (e.g.Smith et al. 2021), but this information is unfortunately not recorded by our main simulations.Instead, we save the gas density and galactic spherical radius at which stellar particles are born, providing a useful proxy for the distribution of young stars that are the main feedback contributors.We show their distributions in Figure A1.
First comparing runs with an increasingly resolved cooling length (top panels), we observe significant differences between their gas density and radii distributions of stellar particles at birth.We compute the two-dimensional Kolmogorov-Smirnov (KS) distance from the reference distribution ('ISM 18 pc') and find (0.06, 0.10, 0.08) for the density distributions and (0.09, 0.10, 0.09) for the galactic radii.Such statistics would be enough to reject that the two samples are drawn from the same underlying distribution at high confidence (p-values < 10 −7 ).
This rejection, however, likely arises from the large sample size (≈ 10, 000) affecting the KS test, rather than from intrinsic differences.In fact, the level of variations between distributions with different cooling lengths is comparable or smaller to run-to-run scatter at a fixed numerical setup (middle and bottom rows).In these cases, we respectively find KS distances of (0.07, 0.09) and (0.04, 0.06) between distribution of gas densities, and (0.06, 0.04) and (0.06, 0.10) between distributions of radii.We also verified that the Wasserstein distances between distributions, which are less sensitive to extrema values, are comparable between stochastic re-simulations and different resolution runs.
We thus conclude that introducing additional refinement on the cooling length does not impact the star formation within our galaxies further than natural stochasticity, with star formation proceeding at a consistent range of densities and radii between all runs.This could have been expected since we minimally modify the AMR structure of the dense disc where star formation occurs.Even when additional resolution is visible in the galactic plane close to the ISM resolution (Figure 2, right-most), such gas is diffuse and thermally-unstable rather than dense and star-forming, and is being better resolved due to its small cooling length.Our star formation criterion based on convergent flows ensures that star formation proceeds in a similar range of densities and radii in this case (green) than in the other runs.To provide an additional verification that feedback conditions in the ISM are comparable between runs, we store the densities within which SNe explode for our four limited-time a-posteriori resimulations described in Section 5.1.Figure A2 shows their distribution across the four resolutions for the cooling length.Differences in these distributions are minimal, giving strong confidence that the star formation, feedback and outflow launching conditions within the ISM are comparable at all resolutions of the cooling length.Rather than differences in initial energetics, the results in Section 3 stem from better resolving the propagation of the outflow into the CGM.

APPENDIX B: SENSITIVITY OF OUR RESULTS TO STOCHASTICITY
Although all simulations in the main suite form a similar amount of stars overall (Figure 3), their instantaneous star formation rate and their ability to drive powerful outflows at a given time also depends on more stochastic parameters, such as the instantaneous ISM conditions and clustering of SNe.In this Appendix, we perform additional simulations to understand the level of variance from different star formation realizations, and confirm that results in Section 3 are robust to this stochasticity.
To this end, we first simulate the base fiducial setup (i.e.only with the quasi-Lagrangian refinement scheme) three times, changing only the random seed of star formation to obtain a different evolution at fixed physical setup.We then repeat this exercise with refinement on the cooling length down to 36 pc. Figure B1 then compares their distribution of star formation rates (top panels), mass and energy outflow rates at 1 kpc (middle rows), and mass and energy outflow rates at 5 kpc (bottom rows) as defined in Section 3.3.As in Section 3.3, we remove inflows from the presented statistics which only occur for the fiducial setups, making trends when comparing to higher resolution runs more conservative.
Focussing first on each suite of three re-simulations (middle and right columns), we observe visible shifts in the median star formation, mass and energy outflow rates that confirm the presence of stochasticity at fixed numerical setup.Nonetheless, star formation rates have well overlapping 16-84 confidence intervals (black lines) both within each individual suite (top middle or top right), across the two stochastic re-simulation suites (top middle and right), and across the resolution study (top left).This further supports that star formation is regulated similarly across all simulations (see also Appendix A).
Mass and energy outflow rates close to the disc (|| = 1 kpc; second and third row) are also very stable against stochasticity, with both re-simulations suites (middle panels, middle and right columns) displaying strongly overlapping 16-84 confidence intervals that also overlap with the resolution study (middle, left).
Away from the disc ( = 5 kpc, bottom rows), mass and energy outflow rates show larger differences between stochastic re-simulations (yellow and red), whether in the fiducial case or when refining on the cooling length (middle and right columns, respectively).We also note that higher average star formation rates (e.g.top, yellow) does not necessarily translate to more vigorous outflows on average (bottom, yellow), reflecting the complexity of sustaining galactic winds.Despite these larger shifts, the stochasticity probed in our fiducial setup (middle column) is incapable of reproducing the sharp increase in outflow loading factors at 5 kpc observed when increasing numerical resolution out-of-the-plane (left and right columns).Even our most vigorous fiducial case (middle column, red) has mass and energy loading factors a factor ≈ 3 short of the quietest run at 36 pc (right We show the average angle of the velocity with the z-direction weighted by specific energy, for the same gas as in Figure 6.Outflows at the highest resolution (cyan) are systematically more open (dotted line showing medians, see also Figure 4), in both the upper and lower outflows (solid and dashed, respectively).But there is considerable scatter around these medians following the varying outflow conditions (see also Figure 3), making a quantitative and conclusive link difficult to establish.

APPENDIX C: OUTFLOW OPENING ANGLE AS A FUNCTION OF RESOLUTION
Figure 4 shows a greater opening of the velocity streamlines with resolution at a single time.Here, we quantify how outflow opening angles vary as a function of time and respond to increasing resolution.
To illustrate this, we select gas in upper and lower cylinder above and below the galaxy as in Section 3.1 ( ≤ 4.0 kpc and 0.5 ≤ || ≤ 5 kpc).For each gas cell, we then compute the angle of their velocity vector and the z-direction.Figure C1 then shows the time series of the average angle for this gas, weighted by specific energy (to upweight the hot and fast gas).Solid and dashed lines show the upper and lower outflows, and dotted lines show the median at each resolution.
We find that, at the highest resolution (cyan), the velocity structure is systematically more open (median towards larger angles) in line with the qualitative findings of Figure 4 and aligning with the increased outflow energetics (Section 3.2).However, there is considerable scatter in opening angles as a function of time, with all different resolutions producing more or less open velocity structures with varying outflow conditions.In fact, the median opening angle does not respond linearly with improving cooling-length resolution ('Fiducial' in purple has the second-most open outflow, but worse off-the-plane resolution).Furthermore, differences between upper (solid) and lower (dashed) outflows are comparable to the shifts with resolution.A causal link between our improved numerical treatment and systematically more open velocity streamlines thus remains tentative, and we hypothesize that repeating this analysis with a longer time baseline, or in systems with higher mass outflow rates might help establish this result.

APPENDIX D: DIFFERENT OUTFLOW GEOMETRIES AT 5 KPC
Here, we test the importance of our chosen biconical shell geometry at 5 kpc.We first show in Figure D1 the same as in Figure 8, but using the full radial shell rather than restricting on polar angles as for the fiducial case.
Figure D1 recovers all the trends observed in Figure 8, in particular the energetic hot outflowing phase being hotter and faster moving at increasing resolution.But using the full radial shell highlights another feature: slow moving gas spreading between two defined tracks at  = 10 4  and  ≈ 5 × 10 6 .Rather than the outflow originating from the galaxy, this gas is part of the hot hydrostatic halo of our galaxy, as evidenced by its sharp drop-off at  =  circ (we also verified that this gas defines a linear track in the T- plane).Our choice of hot halo conditions are fully dictated by the choice of initial conditions for our isolated galaxy set-up, and its contribution can be significant when the diffuse gas is poorly resolved (bottom, left).In our fiducial analysis, we thus aim to mitigate its importance by selecting a biconical shell to more cleanly separate this gas from SNe-powered outflow which we aim to study.
In Figure D2, we verify that trends of outflow loading factors with resolution are unaffected by this choice.We show the time series of the outflow mass and energy loading factors (top and bottom) at each of our resolutions (left to right) computed through our fiducial 2-kpc thick shell at 5 kpc with and without its biconical definition (orange and red), a thinner 0.8-kpc shell (gold), and a 2-kpc thick shell at 10 kpc (brown).
In all cases, the outflow loading factors climb with improved resolution (left to right), validating that our trends are robust to these definitions.As expected due to geometrical factors, the mass loading factors defined through the biconical shell are always lower than the full shell, and mass loading factors decrease between 5 and 10 kpc (top, orange versus purple).However, perhaps surprisingly, the energy loading factors increase with radius (bottom) despite the lack of off-the-plane energy injection from stars (e.g. from runaway stars).This arises due to the rising contribution from the hot  ≈ 10 6 K hydrostatic halo, that acts as heat bath in which the galaxy is embedded.This hot halo can still be pressurized by the SNe-powered outflows and slowly expand, but most of its properties are dictated by our choice of initial conditions.To avoid this issue, we thus select a biconical geometry for our fiducial analysis in this study, and plan to extend our work to a cosmological context studying better-resolved inflows and outflows in a self-consistent environment.

Figure 1 .
Figure1.Visualization of the impact of better resolving the diffuse gas cooling length (right) compared to a traditional ISM-focussed quasi-Lagrangian strategy (left).Panels show edge-on maps of the gas cooling length (top) and its ratio with the spatial resolution of the simulation (bottom), density-weighted along the line of sight at a time of comparable outflow rates ( = 540 Myr, ≈ 0.05 M ⊙ yr −1 ).When resolution is focussed on the ISM (left), the dense outflowing material (top, orange to white,  cool ≈ 10 − 100 pc) mixes rapidly into a diffuse hotter halo (blue,  cool ≥ 100 pc).As the resolution degrades rapidly in these low densities (e.g.Figure2), cooling processes are marginally or significantly under-resolved (bottom left, blue).Forcing refinement on  cool (right, here down to 18 pc) resolves the thermodynamics of the outflowing gas much better, leaving only thin interfaces of under-resolved cooling gas outside the disc plane (bottom right).This improved treatment significantly enhances outflow energetics (Section 3) and modifies their ionic structure (Section 4).Accurately capturing the launching and cooling of outflows within the dense ISM requires sub-pc resolution (top, dark red) posing a distinct, complementary challenge for galaxy formation simulations.

Figure 3 .
Figure 3. Stellar mass growth as a function of time (top), star formation rates averaged over 10 Myr (middle) and mass outflow rates at |  | = 1 kpc(bottom) for our four galaxies with increasingly resolved outflows.Despite the change in numerical setup, all galaxies form a similar amount of stars, ensuring that the total feedback budget is approximately the same across all simulations.Better resolving the cooling length also yields comparable average star formation rates and mass outflow rates close to the disc, providing a controlled study.We discard the first 200 Myr of each simulation to allow the artificial transient phase induced by our initial conditions to dissipate.

Figure 4 .
Figure 4. Thin slices through the gas density (top) and temperature (bottom) of the galaxy viewed side-on at the same time as Figure 2. Increasingly resolving the cooling length of the diffuse gas (left to right) yields a visibly more structured outflow, better capturing the density and temperature contrasts between the hot, volume-filling cavities after shocks and the thin, dense shells of expanding material.At this time, velocity streamlines (bottom, black arrows) are also increasingly open, indicating a faster moving and more energetic outflow.Geometrical interfaces through which we quantify outflow energetics are visualized in the top panels (white lines).

Figure 5 .
Figure5.Time-averaged, mass-weighted distributions of density (top) and temperature (bottom) of the outflowing gas (0.5 ≤  ≤ 5 kpc; (  ) = ()) across our suite of simulations.A better resolved cooling length (purple to green) leads to longer tails towards both denser and more diffuse gas, and both colder and hotter gas.This points towards an increasingly multiphase structure of galactic outflows with increasing resolution (Figure6).

Figure 6 .
Figure6.Time-averaged, mass-weighted temperature-density (top) and pressure-density (bottom) distributions of the outflowing gas when increasingly resolving the cooling length (left to right).In all cases, most gas mass is warm and in photo-ionization equilibrium with the UV background (≈ 10 4 K).However, an increasing resolution yields a better-sampled diffuse phase with an increased scatter in temperatures and pressures at a given density, and allows denser and colder gas at high densities.This increasingly multiphase structure in turn affects how the thermo-kinematics of outflows as they propagate from close to the disc (Figure7) into the inner CGM (Figure8).

Figure 7 .
Figure 7. Time-averaged temperature-velocity distributions with increasing resolution (left to right), weighted by mass (top) and specific energy (bottom) for gas within a 0.1-kpc thick slab at |  | = 1 kpc.Most outflow mass (top panels) is warm and moving slower than the escape velocity of the system (grey dotted in top panels).But better-resolved outflows are systematically more energetic (bottom), extending towards hotter temperatures and faster moving velocities (top-right corner of each bottom panel) along a track of near-constant Mach number (grey dashed).Each level of increased resolution also makes the cloud of hot and subsonic gas more prevalent, indicating that we increasingly resolve the sonic transition of the hot phase.This is key to accurately capturing the acceleration of the hot gas and its subsequent energetics as it escapes into the CGM (Figure8).

Figure 8 .
Figure 8. Same as Figure7but for gas within a 2-kpc wide shell at  = 5 kpc with opening angle  = 45 • along the minor axis of the galaxy (see Figure4for a visualization).At these outer radii, outflow mass is hotter and faster moving ( ≥ 10 km s −1 ) than at |  | = 1 kpc.The improved numerical treatment of the thermo-kinematics of the gas significantly boosts its energetics, with outflow mass systematically travelling at faster velocities (top panels) and a more energy-loaded hot-fast phase (bottom panels).The combination of these factors translate into large boosts in outflow loading factors (Figure9).

Figure 9 .
Figure 9. Mass, energy and metal loading factors (top, middle and bottom rows, respectively) measured through a slab close to the disc (|  | = 1 kpc,  = 4.0 kpc; left) and a biconical shell at large radii ( = 5 kpc,  ≤ 45 • ; right).Violins show the distribution of loading factors over time as the star formation fluctuateduring each galaxy's evolution (median in diamonds and 16-84 confidence interval as black line).Better resolving the dynamics of multiphase outflows from the low resolution fiducial case (purple) to the highest resolution case (green) systematically increases all loading factors.Close to the disc (left), increases are limited for the mass and metal loading factors, but grow to a two-fold increase at larger radii (right).The better-captured propagation of the hot phase significantly enhances the energy loading factors through both interfaces, culminating in a 5× boost at 5 kpc.

FrequencyFigure A1 .
Figure A1.Distribution of gas densities (left) and galactic radii (right) for star particles at their time of birth.Differences between runs with increasingly resolved cooling lengths (top row) are comparable with run to run scatter at a fixed setup (middle and bottom rows), confirming that star formation proceeds similarly in all simulations.

Figure A2 .
Figure A2.Distribution of gas densities at which SNe explode in the four re-simulation tests described in Section 5.1.Differences in distributions are minimal, confirming that SN feedback proceeds similarly at all resolutions of the cooling length.

Figure C1 .
Figure C1.Time evolution of the opening angle of the hot outflowing gas.We show the average angle of the velocity with the z-direction weighted by specific energy, for the same gas as in Figure6.Outflows at the highest resolution (cyan) are systematically more open (dotted line showing medians, see also Figure4), in both the upper and lower outflows (solid and dashed, respectively).But there is considerable scatter around these medians following the varying outflow conditions (see also Figure3), making a quantitative and conclusive link difficult to establish.