Cloud Atlas: Navigating the Multiphase Landscape of Tempestuous Galactic Winds

Galaxies comprise intricate networks of interdependent processes which together govern their evolution. Central among these are the multiplicity of feedback channels, which remain incompletely understood. One outstanding problem is the understanding and modeling of the multiphase nature of galactic winds, which play a crucial role in galaxy formation and evolution. We present the results of three dimensional magnetohydrodynamical tall box interstellar medium patch simulations with clustered supernova driven outflows. Fragmentation of the interstellar medium during superbubble breakout seeds the resulting hot outflow with a population of cool clouds. We focus on analyzing and modeling the origin and properties of these clouds. Their presence induces large scale turbulence, which in turn leads to complex cloud morphologies. Cloud sizes are well described by a power law distribution and mass growth rates can be modelled using turbulent radiative mixing layer theory. Turbulence provides significant pressure support in the clouds, while magnetic fields only play a minor role. We conclude that many of the physical insights and analytic scalings derived from idealized small scale simulations translate well to larger scale, more realistic turbulent magnetized winds, thus paving a path towards their necessary yet challenging inclusion in global-scale galaxy models.


INTRODUCTION
Understanding the inherently nonlinear, dynamical structures that underlie the complex ecosystems of galaxies is a challenging task.It is, however, essential towards making the most of our evergrowing body of increasingly detailed observations.A multitude of open problems still surround galaxy formation and evolution today.One that sits at the very heart is the challenge to understand the multiphase nature of not just the gas within galaxies and their surrounding environment (circumgalactic medium (CGM)), but also the flow of this material in and out of galaxies and how this regulates galaxy evolution (often dubbed the Cosmic Baryon Cycle; see Tumlinson et al. 2017).Galactic outflows driven by feedback mechanisms carry material outwards (Veilleux et al. 2005) while inflowing gas provides fuel for new star formation (e.g.Kereš et al. 2005;Dekel & Birnboim 2006;Fraternali et al. 2015).This cycling connects processes on stellar (∼ pc) scales to galactic (∼ kpc) and intergalactic (IGM; ∼ Mpc) scales, weaving them into a tightly interdependent and highly multiscale tapestry.
Galactic winds, which play a central role in regulating the baryon cycle, have garnered significant attention in recent years due to the clear evidence of their influence on galaxy evolution.They are thought to be primarily driven by feedback processes originating from either massive stars (for lower mass haloes) or active galactic nuclei (AGN; at the higher end of the mass spectrum), and are necessary ingredients for any realistic model of galaxy evolution ★ E-mail: zunyibrent@ucsb.edu(see Somerville & Davé 2015;Naab & Ostriker 2017, for detailed reviews).The result is the formation of large-scale outflows that transport material out of the galaxy while also shaping the CGM (Tumlinson et al. 2017).
Observations of galactic winds reveal them to be ubiquitous across star forming galaxies (Martin 1999;Rubin et al. 2014).They also exhibit a complex, multiphase structure consisting of cold, cool, warm, and hot gas components (see Table 1 for the temperature ranges we use for each in this work), with the various phases spanning a range of different properties and dynamics (Veilleux et al. 2005;Strickland & Heckman 2009;Steidel et al. 2010;Rubin et al. 2014;Heckman & Thompson 2017;Bolatto et al. 2021).A better understanding of the formation, survival, and growth of cool gas clouds which make up the most readily observable phase via emission and absorption lines, as well as their interactions with the other gas phases, is needed in order to provide vital insights into the mass transport mechanisms and overall energetics of these outflows.This in turn is essential for constructing a comprehensive picture of the multiphase nature of galactic winds and their impact on galaxy evolution.
In parallel, theoretical modeling and simulations of galactic winds have become increasingly sophisticated, aiming to reproduce and understand the multiphase nature and complex dynamics observed in these outflows.Chevalier & Clegg (1985) introduced an early analytic model where mass and energy injection into a spherically symmetric region powered a hot outflowing wind.Building on this simple model, recent theoretical work has included radiative cooling and gravity (Thompson et al. 2016), non-spherical expansion (Nguyen & Phase Temperature Range Cold  < 5 × 10 3 K Cool 5 × 10 3 K <  < 2 × 10 4 K Warm 2 × 10 4 K <  < 5 × 10 5 K Hot  > 5 × 10 5 K Table 1.Temperature ranges defined for the four phases of interest referred to throughout this work. Thompson 2021), more realistic non-uniform injection (Bustard et al. 2016;Nguyen et al. 2023), and frameworks for coupled multiphase evolution (Huang et al. 2020;Fielding & Bryan 2022).In particular, the multiphase models describe analytic prescriptions for evolving unresolved clouds in the cool phase and their coupling with the local hot background wind.
In this work, we focus on winds driven by stellar feedback, primarily in the form of energy released by core-collapse supernovae (SNe).While large scale cosmological simulations are now approaching high enough resolutions where multiphase structures (on scales of order 100 pc) are seen to develop self consistently in the CGM (Nelson et al. 2020;Ramesh et al. 2023) and in galactic winds (e.g., in TNG50 Nelson et al. 2019 andFIRE Anglés-Alcázar et al. 2017;Pandya et al. 2021), the detailed properties of multiphase outflows are still most clearly seen in high resolution galaxy scale simulations (e.g., Schneider et al. 2020;Steinwandel et al. 2022;Rey et al. 2023), or interstellar medium (ISM) patch simulations (e.g., TIGRESS Kim et al. 2017;Kim & Ostriker 2018;SILCC Gatto et al. 2017;Rathjen et al. 2023).The ∼ pc-scale resolution of these simulations makes it possible to resolve the energy injection by SNe and the resulting interactions with the surrounding ISM (e.g., de Avillez 2000; Joung & Mac Low 2006;Hill et al. 2012).Such simulations have also revealed a trove of complications surrounding the launching of these winds.For example, the efficacy of this process is sensitive to the spatial distribution of SNe (Creasey et al. 2013;Martizzi et al. 2015;Li et al. 2017;Smith et al. 2021) as well as spatiotemporal clustering (Kim & Ostriker 2017;Fielding et al. 2018) and self-consistency of treatments of star formation and feedback (Kim et al. 2017;Kim & Ostriker 2018).
The survival and growth of cool gas clouds entrained in these winds has also been a hotbed of analytic and numerical studies in recent years.Theoretically, this picture was initially problematic as these cool clouds should be quickly destroyed by hydrodynamic instabilities during this process (Klein et al. 1994).The timescale for the acceleration of an adiabatic cool cloud is a factor of ≳10 times longer than the timescale for its destruction via Kelvin-Helmholtz (KH) and Rayleigh-Taylor (RT) instabilities (Zhang et al. 2017).This was also verified in simulations (Cooper et al. 2009;Scannapieco & Brüggen 2015;Schneider & Robertson 2017).However, it was recently shown that under the right conditions, these cool clouds can not only survive but even grow within the hot wind (Marinacci et al. 2010;Armillotta et al. 2017;Gronke & Oh 2018, 2020).This is possible when the cooling time of mixed gas is shorter than the destruction timescale.Subsequent studies have shown that exact characterization of this parameter space is complex (Li et al. 2020;Sparre et al. 2020;Kanjilal et al. 2021;Farber & Gronke 2022;Abruzzo et al. 2022a).The mechanism that drives this growth is the formation of a long tail structure of cool gas where mixing and cooling takes place efficiently in turbulent radiative mixing layers at the interface of the two gas phases (Ji et al. 2019a;Fielding et al. 2020;Tan et al. 2021).Crucially, this requires long simulation domains that were not captured in earlier simulations, highlighting the importance of capturing the appropriate scales and boundary conditions that define such problems.The interplay between the cool and hot gas phases and the role of various physical processes are hence crucial in determining the fate of cool gas clouds and their impact on the overall dynamics of galactic winds (Fielding & Bryan 2022).Aside from radiative cooling, other important processes at play include turbulence (Schneider & Robertson 2017;Banda-Barragán et al. 2019), magnetic fields (Gregori et al. 1999;McCourt et al. 2015;Grønnow et al. 2018;Hidalgo-Pineda et al. 2023), and thermal conduction (Kooĳ et al. 2021;Brüggen et al. 2023).Despite these advances, many aspects of galactic wind simulations still rely on simplified models and assumptions, such as idealized cloud geometries or time-constant wind properties.The next step forward is to study the problem under more realistic conditions, incorporating the diverse range of physical processes and interactions that govern the evolution of galactic winds and their multiphase structures.
To that end, in this work we connect results from small scale idealized cloud crushing simulations to large scale multiphase galactic winds by studying how clouds form, evolve and interact with the wind.This mesoscale approach bridges the microscale phenomena controlling individual clouds to the macroscale phenomena controlling galactic scale evolution by simultaneously resolving the lifespan of thousands of individual clouds as well as the driving of the highly turbulent, temporally variable, hot wind that they interact with.The basis of our study builds upon a series of simulations following a similar design to those presented in Fielding et al. (2018), which capture the dynamics in a vertically stratified (512 pc) 2 patch of a galactic disc.In these simulations, clustered SNe lead to the formation of superbubbles which propagate quickly enough to break out of the galactic disc, providing a pathway for the energy and momentum released to vent into the halo in the form of a galactic wind.The shredding and entrainment of cool dense clumps (which span a wide range of scales) is seen in these winds, and is expected to be important in influencing the structure and dynamics of the outflow.We extend the box asymmetrically so as to follow the clouds in the outflow out to a larger height, as well as include additional physics such as self-developed magnetic fields in the turbulent ISM. Figure 1 shows a rendering from our simulation of cool structures (in blue) in the hot outflowing wind (in orange).Our main focus is on analyzing the clouds in these winds.More specifically, we investigate how they are formed during superbubble breakout and their properties relative to predictions from idealized small scale simulations of turbulent mixing layers and individual clouds.
The outline of this paper is as follows.In Section 2, we provide a detailed description of our simulation setup.In Sections 3 and 4, we present simulation results, with the former focusing on bulk properties of the outflow and the the latter on analysis of cloud properties including size distribution and growth.Lastly, we discuss some implications and complications in Section 5 before summarizing and concluding in Section 6.

METHODS
In this work, we carry out 3D magnetohydrodynamical (MHD) simulations using the publicly available code Athena++ (Stone et al. 2020) aimed at capturing the multiphase dynamics in a vertically stratified patch of a galactic disc.SNe corresponding to a star cluster are seeded by hand at the center of the midplane of the disc, which leads to the formation of a superbubble that propagates rapidly enough to break out of the disc and release energy and momentum into the halo in the form of a galactic wind.This builds upon the design of a similar suite of simulations presented in Fielding et al. (2018) where additional details and discussion can be found.In our analyses in the following sections, we focus on studying the formation and dynamics of a resulting population of cool material that gets sown into the wind.Before that, we describe in this section our simulation setup, implementations of various physics, and the initial conditions we adopt.

Setup
All simulations are run in three dimensions on Cartesian grids using the HLLD approximate Riemann solver and the third-order accurate Runge-Kutta time integrator.Our simulation setup consists of a rectangular box with dimensions 512 × 512 × 2048 pc and a fiducial resolution of 2 pc on a uniform grid (we also discuss in this work a simulation run with a higher resolution of 1 pc but which fails to drive a successful outflow).The box is asymmetric-the mid-plane is located a quarter box length from the bottom-so as to follow the outflow on the upper side out to a greater height.We use periodic boundary conditions on the sides of the box, along the  and  directions.We adopt outflowing boundary conditions (zero gradient with inflow explicitly disallowed) along the  direction at the top and bottom.We also implement a density floor of  = 10 −8 cm −3 .

Source Terms
We include a background gravitational profile, driven turbulence, optically thin radiative cooling, photoelectric heating, and SNe injection.Here, we provide more details about the implementation of each of these in turn.We do not include any explicit viscosity or thermal conduction.A  = 5/3 equation of state is adopted throughout.

Gravity
We first consider the density profile of a thin isothermal disc in hydrostatic equilibrium.We begin by assuming that the self-gravity of the gas is a subdominant component of the total gravitational potential and is hence unimportant (i.e., ∼ 2Σ gas  is small, which is true for Σ gas ≪ 1000 ⊙ /pc 2 ).We hence neglect self gravity in our simulations.The spherical potential at a vertical height  above the disc and a radial distance  from the center is and hence the gravitational acceleration is For a thin disc,  << , where Ω =  cir / is the angular velocity.Our equation of hydrostatic equilibrium is thus With  =  2   (  is the isothermal sound speed, related to the adiabatic sound speed by  2  =  2  ), we can easily solve this for the density profile: (5) where  ≡   /Ω is the scale height and  0 is the mid-plane density at  = 0. We can however drop the thin disc assumption and solve for the profile out to larger .We now have Solving this gives us: These are the profiles we adopt in our simulations.We test systems in hydrostatic equilibrium with this gravitational profile in Appendix A.
We note that there are caveats to the simplified approach we have taken of assuming a simple external spherical potential which would affect the profiles above.A more sophisticated treatment should account for both disc and stellar spheroid components (e.g., Strickland & Stevens 2000), turbulent support (e.g., Cooper et al. 2007), a dark matter halo component (e.g., Li et al. 2017) and self-gravity.These are particularly important if one is going to distances ≫ 1 kpc (Li et al. 2017).However, our simplified assumption is sufficient for our box size.

Stratified Turbulence
We drive mixed turbulence on large scales so that the mass-weighted velocity dispersion  ≡ ⟨ 2 ⟩ 1/2 ≈ 10 km/s, which is roughly the sound speed of 10 4 K gas and consistent with observed velocity dispersions in the ISM.Physically, this turbulence is driven by either star formation feedback or gravitational instability in ISM gas.The turbulent kinetic energy injection rate is thus where  2 box is the horizontal box length.The turbulence is driven on the scale of the disc scale height with power distributed evenly between wave numbers  = 2 to  = 3.The turbulence is constrained to follow a Gaussian profile with scale height  ∼ 80 pc (we test this implementation of constrained turbulence in Appendix B).We drive the turbulence every 5 × 10 −3 crossing times  cross .The driven turbulence in Athena++ adopts an Ornstein-Uhlenbeck process to smoothly vary the velocity perturbations over some correlation time ( Lynn et al. 2012), which we set to 8 Myr.Turbulence along with heating and cooling leads to the formation of a turbulent multiphase ISM which we allow to form for 60 Myrs (roughly a turnover time) prior to any supernova explosions.

Cooling/Heating
The net cooling rate per unit volume is formulated as where Λ is the cooling function and Γ is the heating rate.Our cooling function Λ() combines the collisionally ionized equilibrium (CIE) cooling curve in Gnat & Sternberg (2007) for  ≥ 10 4 K with the cooling function for  ≤ 10 4 K in Koyama & Inutsuka (2002).We obtain our cooling curve by performing a piece-wise power law fit over ∼ 50 logarithmically spaced temperature bins, starting from a temperature floor of 10 K up to a maximum temperature of 10 10 K. We also include a photoelectric heating (PEH) rate Γ = 10 −26 (⟨  ⟩/ cm −3 ) erg s −1 , where we have scaled the solar value by the average mid-plane density to approximate the scaling of PEH with higher star formation rates in denser regions.This approach is similar to those used in Kim & Ostriker (2017) and Fielding et al. (2018), which run similar setups.We assume solar metallicity ( = 0.7,  = 0.02).We also use a fixed mean molecular mass  ∼ 0.62 for a fully ionized plasma.This has a ≤ 2 factor of error in the temperature of neutral/partially ionized gas below 10 4 K, but should should not affect our overall conclusions.Figure 2 shows the aforementioned cooling curves along with our fit and heating curves for ⟨  ⟩ = 100 cm −3 at several different pressures (e.g.,  4 is  = 10 4  B cm −3 K).We also adopt an additional constraint on the simulation timestep over the standard CFL constraint where we require that that the timestep is less than or equal to one quarter of the shortest singlepoint cooling time  cool across the whole domain.This is to ensure good coupling between cooling and the hydrodynamical evolution.
Finally, we develop and implement an extended version of the fast and robust exact cooling algorithm described in Townsend (2009) to include heating (See Appendix C for full details).

Supernova Injection
We include both core-collapse supernova and Ia rates using the piecewise power law fits given in Hopkins et al. (2023).The combined SNe rate (shown in Figure 3) is where 2 ,  3 ) = (0.39, 0.51, 0.18, 0.0083) and ( 1 ,  2 ,  3 ) = (3.7,7.0, 44) Myr.The time of the first Ia is set to be after the time of the last core-collapse SNe ( =  3 ).For  ≡  * / ⊙ , the SNe number per initial solar mass  SN ≡  SN / is thus given by The time where  supernovae have gone off is hence We deposit 10 51 ergs per SN as pure thermal bombs.While one could determine the amount of thermal and kinetic energy injected by each SN (see Martizzi et al. (2015), who calibrated their model to high resolution simulations of individual SN remnants), this would matter only for the first couple of SNe.Since the SNe are clustered tightly both in space and time, most of them barring the first few occur in the hot and low density remnant of the previous SNs, assuming that Δ SN <  PE ( PE being the timescale over which the SNR reaches pressure equilibrium or mixes with the ambient ISM) so that a coherent bubble can be driven (Fielding et al. 2018).Hence their cooling radii  cool are around an order of magnitude larger than the injection radius  inj , which means that most of the energy from the SN is transferred to the surrounding gas.
We assume a spherical geometry for the injection site with radius  inj = 20 pc.We assume a stellar ejecta mass of 8.72 ⊙ for corecollapse SN and 1.4 ⊙ for Ia.We also assume that all SN go off at the origin, since the cluster radius  cl = 10 pc is relatively small.
While out of the scope of this work, we note that in reality, stars far from the disc center could have a disproportionate effect in driving winds.For instance, Ia at late times from stars that settle high above the disc or OB (binary) runaways (e.g., Steinwandel et al. 2023).This is mainly due to the dual effects of a larger  cool in the lower density region along with the lower scale height contributing to a large  cool / ratio.
For each cell with center within the sphere radius, we then inject energy from the SN uniformly.Note that there will be some error due to resolving the volume of the sphere, but the error is small as long as we resolve the injection radius by more than a few cells.

Initial Conditions
We adopt a circular velocity  cir = 175 km/s and galactic radial distance  = 1 kpc, which corresponds to a scale height  = 66 pc.
The box is initially filled with  = 10 4 K gas with a mid-plane density  = 100 cm −3 at  = 0.The average density within one scale height is ⟨⟩ = 86 cm −3 .Using ⟨⟩ ≡ Σ g /(2  p ), this corresponds to a gas surface density Σ g ≈ 175  ⊙ pc −2 .For a a star formation efficiency  * ≡  cl /ℎ 2 Σ  ∼ 0.1 (as found in Grudić et al. (2018)), this corresponds to a star cluster mass  cl ∼ 2.5 × 10 5  ⊙ .For comparison, the simulation from the suite of TIGRESS simulations under the SMAUG project with the highest star formation rate (R2) has Σ g ≈ 74  ⊙ pc −2 (Kim et al. 2020a).We initialize a magnetic field aligned along the x-direction with plasma beta  ≡  gas / B = 100 everywhere.However, through cooling and turbulent dynamo amplification,  ∼ 1 is rapidly achieved in the ISM.

RESULTS : BREAKOUT AND WIND OUTFLOW
In this section, we provide context for the winds that contain our clouds, setting the stage for our main focus-a closer look at cloud properties-in the next section.There is already a large body of existing work studying these SN driven outflows-we concentrate mainly on analysis meant to highlight the most salient points that are relevant to these clouds.

Superbubble Breakout
Supernovae (SNe) are powerful energetic events that release a tremendous amount of energy (∼ 10 51 ergs) and mass into their surrounding environments.However, individual SNe are ineffective at driving galactic winds because most of this energy is radiated away before it can break out of the disc.Spatiotemporal clustering of SNe resolves this problem, and is motivated by the fact that stars form primarily in clusters.When multiple supernova explosions occur in the same region of space and over a relatively short period of time, they can overlap and form superbubbles.These superbubbles are able to drive the expansion of the material outwards from the cluster, breaking out of the disc before the cluster runs out of SNe.After superbubble breakout the hot material inside is able to vent its energy into the halo and power an outflowing galactic wind.As the superbubble propagates through the clumpy multiphase environment, it sweeps around higher-density regions, resulting in a surface that is fractal in nature.This complex interplay leads to the fragmentation and break up of the ISM, seeding the hot wind with a population of cool clouds.Figure 4 illustrates this by showing temperature slices at various stages of the process described above.

Asymmetric Outflows
By virtue of having to make its way through a multiphase ISM, the expansion and eventual break out of the superbubble can exhibit great asymmetry in terms of driving outflows above and below the disc.This can be self-reinforcing-whichever side breaks through first is then free to vent into the low density halo, creating a channel of lower resistance for the hot gas to funnel through and leading to a weak or non-existent outflow on the opposite side.In addition, turbulent motions in the ISM can exacerbate or reverse this asymmetry.For example, if cold dense material in the ISM moves over the cluster region, it 'caps' and inhibits the path of the hot expanding wind, much like turning off a valve, resulting in a redirection towards the other side of the disc.In this manner, an outflowing wind can be quenched at the base on one side midway during the lifetime of the cluster SNe.While we did not investigate this in further detail, we observed    that this was a common occurrence in our simulations-due to the nature of the setup, a wind that breaks out on one side pushes the ISM radially outwards in the plane of the disc.Because of the periodic boundary conditions, this leads to a build up of ISM material on the sides which then falls back towards the cluster.This sloshing motion often eventually halts the wind and drives it to the opposing side.
In a more realistic larger scale setup, this indicates that dynamical interactions between different star clusters could complicate wind driving, either reinforcing or disrupting outflows depending on their spatiotemporal separations.In Appendix E, we show and analyze a case where the wind is quenched.
In Figure 5, we show the outwards mass flux per unit area at four different heights in our box as a function of time for our fiducial simulation.The blue and orange lines show the mass fluxes below and above the disc respectively at 0.5 kpc, while the red and teal lines are at 1 and 1.5 kpc above the disc.The dashed orange line shows when the mass flux becomes negative, i.e., material is falling back to the disc.We can see that the breakout is initially on the top side which launches an outflowing wind.As the outflow is interrupted at the base, the wind loses its power source, and begins to slow and dissipate.Meanwhile, the outflow is directed to the bottom side of the disc.This eventually reverses again nearing the end of the cluster SN period.The red region highlighted in Figure 5 shows the time range in which we analyze clouds in Section 4. This choice corresponds to a time period right after the peak of the SNe rate (see Figure 3).To illustrate this more clearly, Figure 6 shows temperature slices at various times in the same simulation.We see that the wind initially breaks out on the topside and that fragments of the ISM become entrained in the hot outflowing wind, whereas it fails to do the same on the bottom side, which is apparent when looking at the rightmost two panels at later times.

Turbulent Winds
Figures 7, 8, and 9 show the range of quantities we track in our simulations, and allow us to make some general statements about the outflowing winds.Figure 7 shows slices of temperature, outflow velocity, number density, pressure and magnetic plasma beta for a single time snapshot 14 Myr after the first SN, in the middle of the red shaded region shown in Figure 5.This shows the general picture of clustered SNe driving a hot outflowing wind which contains a population of cool clouds, which we focus our analysis on in the Figure 9. Slices of the passive scalars we have used in this simulation.We use a total of 5 passive scalars.The first two are introduced right before the first SN explosion.All cold gas ( < 5 × 10 3 K) is dyed with passive scalar  0 , and all cool gas (5 × 10 3 K <  < 2 × 10 4 K) is dyed with passive scalar  1 .Passive scalars  2 and  3 track how much of the gas in a cell was ever hot or cold.Lastly, passive scalar  4 tracks all SN injected gas.Beyond the disc and the injection site,  2 and  3 are close to unity everywhere, indicating that there is strong mixing between phases.next section.Figure 10 is composed of 2D histograms that show how these properties vary with height.The histograms are volume weighted and hence largely reflect wind properties, especially at heights above 0.5 kpc.The top-most panel shows that the wind has a temperature of  ∼ 10 7 K, along with cool clouds that have a temperature of  ∼ 10 4 K.In the second panel from the top, it can be seen that pressure drops off rapidly away from the disc.The pressure in the wind is roughly /  ∼ 10 4.5 cm −3 K .In terms of non-thermal components in the wind, both magnetic fields and turbulence are small compared to thermal pressure.We can see from the third panel that although the disc is strongly magnetized with plasma beta  ∼ 1, the wind itself has a very weak magnetic field with  ∼ 10 2 − 10 3 .Lastly, the bottom-most panel shows that above 0.5 kpc, the turbulent Mach number in the wind  turb ∼ 0.25.This brings us to a discussion of the nature of turbulence and mixing in the wind, where we will also define how  turb is computed.
Dispersed multiphase flow systems such as particles or droplets in liquid or gaseous flows are common in various engineering fields (Balachandar & Eaton 2010).When the mass fraction of the dispersed phase is comparable to that of the carrier phase, the back reaction is non-negligible and the system is said to have two-way coupling.The nature of this multiphase turbulence is still considered an open problem due to the involved complexities.For example, it is unclear how particles can affect turbulence in the carrier phase-studies have shown generally that small particles can attenuate turbulence by dissipating energy for turbulent eddies, while larger particles can amplify turbulence through mechanisms such as vortex shedding (Gai et al. 2020).Similarly, it is apparent in our simulations that the presence of cloud fragments induces turbulence on the wind scale.It has also been widely shown that instabilities on cloud scales such as KH, RT, and Richtmyer-Meshkov not only disrupt clouds but mass load the wind with turbulent gas (Klein et al. 1994;Mac Low et al. 1994;Pittard & Parkin 2016;Banda-Barragán et al. 2020).This arises firstly during the breakout phase as the superbubbles propagates through the multiphase ISM.Once this reaches the outflow phase, the flow of the hot wind around the population of clouds similarly induces a large scale turbulent field on the scale of the disc scale height.Figure 8 shows slices of the velocities in the three axisaligned directions.The vertical velocities  z reflect the bulk outflow velocity, while the velocities  x and  y parallel to the disc plane reflect the turbulence in the wind.
To examine the turbulence more closely, we apply a Gaussian filter to estimate the turbulent velocity  turb .This method involves approximately removing large scale variations related to the bulk flow.This is done by convolving each velocity component with a Gaussian kernel to estimate the bulk velocity along that component.Subtracting this thus leaves the turbulent velocity of that component.While exact results depend on the kernel size, qualitative conclusions remain unchanged (Abruzzo et al. 2022a).We choose to use 16 pc, which sits comfortably within the range of sizes spanned by the cool clouds.The -th component of the turbulent velocity is hence defined here as where   () is the three-dimensional Gaussian filter.We weight the filter by density so that the velocities within the cloud are not dominated by the that of the hot gas.The right most panel of Figure 8 shows how the turbulent Mach number using the above definition for turbulent velocity ( turb ≡  turb / s ) varies in the slices.In general, the turbulent Mach number in the wind is high only close to the disc, and quickly weakens away from it.On the other hand, the cool clouds have transonic internal turbulence.This is shown more quantitatively in Figure 11, which shows a 2D histogram of  turb and temperature  at the same time snapshot as Figure 10, but only considering heights above 0.5 kpc.The orange dashed line shows the average  turb (∼ 1) within the orange shaded region which represents the cool clouds, while the red dashed line similarly shows the average  turb (∼ 0.25) in the red shaded region which represents the hot wind.In all, there is significant turbulence both within the clouds and in the outflowing wind.This means that these environments are more similar to the turbulent boxes studied in (Gronke et al. 2022) rather than idealized laminar wind tunnel setups.This is critical because it is turbulence that drives mixing between the phases and hence their coupling.Finally, Figure 9 shows slices of the passive scalars we have used in this simulation, and also illustrates the amount of mixing between phases happening in these winds.We use a total of 5 passive scalars.These passive scalars track gas with certain properties.The first two are introduced right before the first SN explosion.All cold gas ( < 5 × 10 3 K) is dyed with passive scalar  0 , and all cool gas (5 × 10 3 K <  < 2 × 10 4 K) is dyed with passive scalar  1 .This is to track how much cloud material was originally cold and cool ISM material.Passive scalars  2 and  3 track how much of the gas in a cell was ever hot or cold respectively, similarly starting from right before the first SN, so as to quantify the amount of mixing.Lastly, passive scalar  4 tracks all SN injected gas, which similar to  0 and  1 , tracks how much cloud material was from the SN.In the following section, we use  0 ,  1 , and  4 to quantify the original temperatures of cloud material.Here, we will just point out how beyond the disc and the base of the jet, both  2 and  3 are close to unity everywhere.This means that all gas not close to the injection point or part of the ISM has at some point been cold or hot.

Mass and Energy Flux Phase Distribution and Cloud Entrainment
Finally, we can look at the fraction of the wind's mass and energy flux carried by the different phases as defined in Table 1. Figure 12 shows 2D slices of the general location or each of these phases at 14 Myr.Cold gas is only located within the ISM.Cool gas is found in the ISM and the clouds that are contained in the wind, which is made up of the volume filling hot phase beyond the disc.Warm gas is found in the interfacial mixing layers.
Figure 13 shows the temperature velocity distributions at different heights above the disc over the time period identified in Figure 5, weighted by mass flux in the top row and energy flux in the bottom row.The PDFs  M (, ) and  E (, ) are weighted by mass and energy flux respectively, where  ≡ log( out ) and  ≡ log() in units of km/s and K.The PDFs  q for each quantity  is defined as 1 ⟨⟩  2   where ⟨⟩ is the average over both time and the horizontal slice, such that ∫  q  = 1.We define the mass flux to be  z and the energy flux to be  z  z where  z = ( 1 2  2 z +  th ) and  th is the specific internal energy.The horizontal grey dashed lines demarcate the boundaries between the phases defined above.Consistent with previous work, the hot phase dominates the energy loading while there is significant mass loading by the cool phase (Fielding et al. 2018;Kim et al. 2020a,b;Rathjen et al. 2023).Mass and energy loading factors are similar to those found in more global simulations such as Schneider et al. (2020), who also observe cold phase mass loading in a volume filling hot wind.This is in line with previous studies of multiphase flows such as the more idealized setups of Banda-Barragán et al. ( 2021), where the hot phase occupies ≥ 90% of the volume but only contains ∼ 10% of the mass.The grey dotted line shows the sound speed of the gas at each temperature-the outflow velocity of the wind is transonic with a Mach number close to unity.We also see that the cool material gets more entrained at increasing heights, indicating that there is momentum transport from the hot to the cool phase, which has been found to be driven by mixing rather than ram pressure acceleration (Melso et al. 2019;Schneider et al. 2020;Tonnesen & Bryan 2021).

RESULTS : CLOUDS
In this section we focus our attention on analysis of the clouds embedded within the hot wind and their properties.We first discuss how we define and identify clouds, followed by how we measure related cloud properties.We then analyze size distributions and various scaling relations associated with cloud growth rates.Finally, we explore other interesting aspects of cloud properties relating to non-thermal support, material origins, and wind alignment.

Cloud Identification
In order to map out and catalogue the properties of clouds embedded in the hot outflowing wind, we must first define them.We identify clouds by applying a threshold at  = 10 5.5 K, represented by the red dashed line the upper panel of Figure 14.Any gas below this threshold is thus considered 'cloud material'.The lower center panel of Figure 14 shows this cool material identified from the corresponding temperature slice.While we use a fixed threshold, we find that cloud identification in our simulations is generally insensitive to the exact value used. 2 We can understand this physically-most cloud gas is at a temperature ∼  = 10 4 K, while most of the wind is at  > 10 6 K. Gas at intermediate temperatures mainly exists only at the thin interfaces between these two phases in turbulent radiative mixing layers.This can be seen visually in the temperature slice in the lower left panel of Figure 14 (also see the 'warm' panel in Figure 12).It is hence not surprising that the identification of clouds is understandably robust to the choice of threshold provided that this threshold temperature lies in the mixing region.The simple approach of using a fixed threshold is thus as effective as it is because of the strong biphasic nature of the wind and the clouds.This can be seen in the upper panel of Figure 14 which shows the temperature distribution in the wind at a single time snapshot, along with the threshold value used.Using a constant threshold also allows for more consistency over multiple time slices.We can thus define a single cloud as a collection of cells that fall below the threshold and that are interconnected (where connected is defined as sharing at least one corner).Except for the analysis on cloud size distributions, where we account for clouds that intersect the periodic dimensions and hence wrap around, we ignore clouds that are touching the borders of our simulation.This means we naturally exclude disc material.We also only consider clouds that have a volume greater than 16 cells.The bottom right panel of Figure 14 shows an example of these clouds arbitrarily colored by their cloud identification number.We also compile our catalog over multiple time snapshots (a total of 10 spanning 8 Myr).

Measuring Cloud Properties
For each identified cloud, we measure a number of properties as follows.The position of the cloud is computed as the volume weighted centroid of cloud material.Individual cloud properties such as velocity v cl , density  cl , thermal and magnetic pressures  cl and  B,cl are volume averaged over cool (< 1.01 × 10 4 K) gas in the cloud.Temperature is density weighted, so that  cl =  cl  cl .Passive scalars are also density weighted.We define the local background environment of a cloud as any hot (> 10 6 K) wind gas within a cube centered on the cloud that has a length on each side equal to twice that of the longest axis of the tight bounding box that contains the cloud, and also require that there are at least 16 such cells in this volume.Corresponding wind properties such as the background wind velocity v w are averaged in the same way as cloud properties.The turbulence in the wind v turb,w is taken to be the rms velocity in the wind frame in this environment.
We adopt a geometric approach towards characterizing the turbulence on the surface of a cloud, which involves measuring quantities on a temperature isosurface constructed using the marching cubes algorithm described in Lewiner et al. (2003).This isosurface is represented by a mesh of triangular faces.Figure 15 shows examples of this mesh in the lower panel for various clouds.The cloud in the upper left is shown in greater detail in the upper panel.
To measure the turbulence in the mixing layer on this surface, we use the method described in Abruzzo et al. (2022a), which they demonstrate provides a measure of turbulence that is comparatively robust to gradients in the laminar component of the background wind as compared to the other approaches of characterizing turbulence that they explored (filtering methods and velocity structure functions).The outline of the method is as follows: Consider a point on the temperature isosurface that is at the centroid of one of the triangular faces which constitute the mesh.Working in the cloud frame, we linearly interpolate the velocity v and the logarithmic density  at that point.We define a set of axes n, ŵ and φ, where n is the inwards direction normal to the surface, ŵ is the direction of the component of the background wind  wind in the surface plane and φ is the direction on the surface orthogonal to both n and ŵ, i.e., φ ≡ n × ŵ.These definitions are physically motivated-the normal direction captures turbulent radiative mixing layer driven accretion, while the wind direction captures the coherent flow of the background wind.Picking out the φ direction thus allows us to disentangle the turbulence from these flows.The velocity component is hence simply v  ≡ v • φ.We thus estimate  turb to be the area weighted standard deviation of v  .The upper panel of Figure 15 shows how v  varies over the surface of a cloud.The top right inset shows a surface area weighted histogram of v  on this surface.As expected, the distribution is normal and centered on zero. turb for this cloud is simply the standard deviation of this distribution.The inflow velocity through this surface is similarly given by v in ≡ v • n.
In Appendix D, we detail how we estimate cloud growth rates , comparing the effectiveness of various methods of doing so.Ultimately, we use the net inwards mass flux on the temperature isosurface (  surface ≡ surface  in Δ) with a scaling factor calibrated to match the total cooling luminosity, so  ∼  scale  surface .Appendix D also presents an important caveat that these methods fail when the cloud is not actually growing, and thus do not reflect the transition to the cloud destruction regime.Instead, we show that we can estimate this scale from the cloud size distribution.Now that we have discussed how we identify clouds and measure their properties, we are ready to further analyze these measured properties.We first present the distribution of cloud sizes in the simulation and compare to analytic formulations.We then compare estimated cloud growth rates with predicted scalings from models for the surface area and associated mass flux as given in equations ( 13) -(20).

Distribution of Cloud Sizes
The snapshots shown in Figures 4 and 6 present an intriguing picture of an initial population of cold clouds seeded in the wind via direct fragmentation of ISM material.Here, we look at the implications of this process on the cloud size distribution.
How are these clouds generated?Unlike setups where cloudlets are formed via gravitational collapse or thermal instabilities, the process of superbubble breakout here itself leads to fragmentation of the ISM.While a detailed study would be needed to fully understand this fragmentation process, we believe there are three key pieces at play.Firstly, the bubble expands through a turbulent multiphase ISM, which was itself generated through a combination of driven turbulence and thermal instability.Both the turbulent field and the multiphase landscape distort the bubble geometry on a large range of scales and can lead to the formation of disconnected cloudlets.In this case, the nature of the turbulence and the initial distribution of cold ISM material might play important roles.Secondly, the thin shells of the expanding bubble interfaces are susceptible to a host of instabilities such as Kelvin-Helmholtz, Rayleigh-Taylor and Richtmyer-Meshkov, and have been shown to themselves be fractal in nature (Fielding et al. 2020;Lancaster et al. 2021).Lastly, the expansion of the bubble through this multiphase inhomogenous medium leads to large scale turbulence (independent of the initial driven turbulence in the ISM).This can also drive turbulent fragmentation (e.g.Kolmogorov-Hinze theory for turbulent fragmentation of droplets; Kolmogorov 1949;Hinze 1955).The detailed influence of each of the above, along with their interplay, will be an important focus for future work.Regardless, the simulations suggest that it is the fragmentation of the initial ISM that determines the distribution of cloud sizes, at least when the outflowing wind first breaks out.The initial cloud mass population in the hot wind can thus be viewed as being seeded by the fragmentation of the ISM, driven by the process of superbubble breakout.We will show below that our simulations support this picture of cloud genesis via ISM fragmentation during superbubble breakout, and leads to a distribution that scales as: where  () is the number of clouds with mass .Inverse power law distributions with a -2 exponent, which indicates that there is a constant contribution from each logarithmic bin, are commonly observed.We refer the reader to Section 5 for further discussion.What sets the lower and upper limits on cloud sizes?Since the clouds are formed during the breakup of the ISM, the largest cloud size is simply constrained by the initial scale of the ISM, i.e., the scale height  of the disc.The scale of the smallest clouds is a more complex question.While fragmentation can lead to arbitrarily small clouds, they should be rapidly destroyed in the hot turbulent wind below some critical size, where growth due to cooling becomes ineffective.This survival threshold  crit presents a natural choice for the lower limit.However, the problem of what determines  crit is a thorny one, and there has been no shortage of recent investigations into the matter.We will briefly review some of them.
For a cloud of size  cl with overdensity  relative to its background and temperature  cl at rest embedded in a hot wind of temperature  hot moving at a velocity  wind , Gronke & Oh (2018) identified a characterize size  crit,cc above which the clouds survive till entrainment and grow.This size corresponds to where the cloud crushing time  cc ∼ √  cl / wind , the timescale on which a cloud is destroyed by hydrodynamic instabilities (Klein et al. 1994), becomes longer than the cooling time of mixed gas  cool,mix ≡  cool ( mix ), where  mix ∼ √  cl  hot (Begelman & Fabian 1990).From their turbulent box simulations, Gronke et al. (2022) found that there is also an additional empirically determined Mach number dependence of the functional form  ( turb ) ∼ 10 0.6 turb .They attribute this additional Mach dependence to increased turbulent disruption, in contrast to the increased survival times seen in cloud crushing simulations with high Mach winds which stems from cloud compression (Scannapieco & Brüggen 2015; Bustard & Gronke 2022).The critical size for survival is thus given by where  cl,4 ≡ ( cl /10 4 K),  3 ≡ /(10 where  hot is the density of the hot gas.At  = 100, the difference between the two criteria is small, and when comparing the two, attributed to different ways of defining cloud survival (Kanjilal et al. 2021).At larger  however, these criteria differs by orders of magnitude.Abruzzo et al. (2022a) find that at  > 100, the Li/Sparre criterion agrees much better with simulation results, but also point out that  cool,wind cannot be physically important since their results are unchanged when they artificially shut off wind cooling.They instead propose a new criterion  cool,mix <  shear which captures the empirically derived  scaling of  crit in a physically motivated model (Abruzzo et al. 2023).Here  shear ∼  cl / wind is the wind crossing time over the cloud.The physical argument is as follows.As hot phase fluid elements travel past the cloud, they mix with the cooler fluid elements from the cloud to produce intermediate temperature gas at the interface.Clouds will only grow if a parcel of gas is able to mix and cool before being advected past the cloud.This requires that the cooling timescale of the mixed gas is shorter than the wind crossing time past the cloud.This criterion gives a cloud survival size Figure 16 shows the cumulative size distribution of clouds at different times in the simulation.The distribution of / ∝  −2 is equivalent to a cumulative distribution that scales as  (> ) ∝  −1 .Instead of a mass scaling, we consider instead the volume scaling  (> )/ since the fragmentation process is not driven by a density dependent process, although the difference made by this choice does not affect the results.The pink region in Figure 16 indicates where clouds are resolved by less than 16 cells (and hence not included in the catalog of cloud properties).The grey region show the −1 power law scaling for comparison.The upper panel shows distributions that start from when the superbubble first breaks out from the disc.The number of clouds in the wind increases over time as the disc continues to fragment.We note that the  (> ) ∝  −1 power law distribution is seen even at early times, supporting the argument that this is a result of the disc fragmentation process during the process of superbubble breakout.We have checked that the the distributions does not show any consistent trends with height in these simulations.
The distribution exhibits a flattening at small scales somewhat above the resolution limit, which is expected for clouds below the critical radii for sustained growth.To estimate these critical radii we take a typical set of parameters for these winds.For example, in our fiducial simulation, these are  wind,rel ∼ 0.4,  turb,wind ∼ 0.5,  3 ∼ 30 and  ∼ 400. Figure 17 shows the distributions of  wind,rel and  along with their means.This gives a cloud survival size of  crit,shear ∼ 5 pc.
This nicely match the scales below which we see the distributions flatten.The vertical grey dashed lines in Figure 16 show the upper and lower bound volume estimates corresponding to  crit,shear and the disc scale height  ∼ 66 pc respectively.Ideally, the smallest surviving clouds should be resolved by multiple cells to achieve first order convergence in terms of total cloud mass.With a resolution of 2 pc, we resolve  crit,shear by a factor of ∼ 2. In Appendix E, we perform the same analysis for a simulation which fails to sustain an outflowing wind.

Theoretical Model
Assuming that the conditions are right for a cloud to survive and grow in a wind, we can write its mass growth rate as where  wind is the density of the hot wind,  in is the inflow velocity corresponding to the mass flux from the hot background onto the cloud, and  cloud is the effective surface area of the cloud.
The first ingredient in this cloud growth model is  cloud .It is important here to distinguish this effective surface area from the non-convergent surface area of the cool gas (Fielding et al. 2020;Abruzzo et al. 2022a), which corresponds more closely to the temperature isosurfaces we measure in our simulations below.Instead, a more accurate characterization of  cloud is an effective surface area corresponding to some smooth envelope around the cloud (Gronke & Oh 2020).In cloud crushing simulations, surface area has an initial rapid growth stage associated with tail formation, followed by slower isotropic growth (Gronke & Oh 2020;Abruzzo et al. 2022a).This corresponds to  cloud going from initially scaling as ∝  (linearly with volume), to scaling as ∝  2/3 .The initial exponential growth can be attributed to a tail formation phase where the high shear mainly leads to mass growth in a growing tail downstream of the cloud.Eventually, the cloud becomes entrained and the surface area scales much as a monolithic cloud would.However, it was found in Gronke & Oh (2023) and Tan et al. (2023) that for clouds which never get entrained (either because of large scale turbulence or gravitational acceleration), and are hence subject to continuous shear, the cloud surface area instead scales as  cloud ∝  5/6 .This scaling can arise from the fractal nature of the mixing surface, which has been measured in mixing layer simulations to have a fractal dimension  ∼ 2.5 (Fielding et al. 2020), which would mean  ∝  /3 ∼  5/6 .If we assume that the smallest surviving clouds are the most spherical, we can model the effective cloud surface area as where  crit and  crit are the spherical area and volumes respectively corresponding to the critical cloud survival size  crit , which we discuss below.
Next we discuss the second ingredient in the model, the inflow velocity  in .In plane parallel simulations of radiative turbulent mixing layers (Tan et al. 2021;Fielding et al. 2020),  in is found to scale as where  turb is the turbulent velocity in the mixing layer,  turb is the largest mixing scale (typically the outer scale of the turbulence), and  cool,min is the minimum cooling time in the simulation, which scales inversely with pressure ( cool,min ∝  −1 ).Additionally, these simulations find that the turbulent velocity scales with the relative shear velocity in these setups, such that  turb ∼  rel  shear , where  rel is some constant of proportionality (Tan et al. (2021)   velocity drives ongoing turbulence after entrainment (Jennings & Li 2021;Abruzzo et al. 2022a).To account for this saturation at low shear velocities, we can hence express the turbulent velocity as In our simulations, we find that  rel ∼ 1/15.The above scalings assume that we are in the subsonic to transonic regime where the relative velocity between the hot wind and the cloud does not exceed the hot gas sound speed.The turbulent mixing velocity saturates as clouds becomes supersonic (Yang & Ji 2023), which would change the scalings above.We also assume that we are in the rapid cooling regime, where the cooling time is much shorter than the turbulent mixing time / turb , i.e., the outer scale eddy turn over time (Damköhler number  ≡ /( turb  cool,mix ) > 1; Tan et al. 2021).
In the following sections, we compare the model outlined above to results from the clouds in our simulations.

Cloud Surface Area
The isosurface method detailed previously allows us to compare the scaling of the cloud surface area as a function of cloud size to the model scaling in equation ( 18).We take the cloud surface area  cl to be the area of the isosurface.We define this to be with respect to a measurement scale of 2 pc when constructing the isosurface, since the isosurface area varies with this choice of scale (Fielding et al. 2020).We do not expect this area to be the same as  cloud from equation ( 18), since its magnitude in part depends on the choice of scale and temperature of the isosurface, but we do expect the scaling with size to remain the same.Figure 18 shows a 2D histogram of  cl and cloud radius  cl (the radius of a sphere with equivalent volume), along with the  cl ∝  5/6 ∝  5/2 cl scaling from equation ( 18).We find that most of the clouds agree with this scaling, except for the largest clouds which appear to have a steeper slope closer to  cl ∝  11/4 cl or even  3 cl .We discuss this further below.

Cloud Growth
We can also see how well the scalings in equations ( 19) and ( 21) hold up.For reference, we expect the inflow velocity  in to scale as in the strong cooling regime where the mixing time  turb / turb is longer than the cooling time of mixed gas  cool,mix (Tan et al. 2021).
For a typical cloud,  mix ∼ 10 5.5 K and  3 ∼ 30, giving  cool,mix ∼ 0.15 Myr.In the mixing layers,  turb ∼  s,cool ∼ 15 km/s.This suggests that we are in the strong cooling regime if  turb > 2 pc.We discuss this in further detail shortly, but  turb is usually taken to be the outer scale of the turbulence, where mixing is the most effective, hence we expect our clouds to comfortably fall within this strong cooling regime.In fact, this criterion from Tan et al. ( 2021) is likely oversimplified since the cooling curve drops off sharply above 10 5 K-it is reasonable to expect that the requirement of  turb might be even smaller.
To compare cloud properties with these scalings, we define following equation ( 17) but using the isosurface area and the associated mass flux through it.Figure 19 demonstrates the  in scaling with,  cl ,  turb , and  w respectively.Here we use  w as proxy for the characteristic mixing layer cooling time  cool since it scales inversely with  w , and we will discuss how  cl relates to  turb .More specifically, each panel show the scaling of  in divided by the expected scalings of the other variables so as to remove any cross-dependencies.Each grey point in the panels represents individual clouds, while red crosses mark a corresponding binned scatter plot to more easily visualize trends in the data (Cattaneo et al. 2019).Binned scatter plots partition the  variable into fewer bins and display the average outcome for observations within that bin.Orange lines show the scalings we expect from our model.The top panel of Figure 19 shows  in /( ) plotted against  cl .In previous work, it has been common to set the outer turbulent length scale  turb to be either the box size (in the case of turbulent mixing layers) or the cloud size (for cloud simulations).In the latter,  turb is in fact set to the initial cloud size and fixed for the rest of the simulation, even if the cloud grows in mass and volume by over an order of magnitude or fragments into many pieces.Here, we expect that if  turb ∼  cl , then we should see that  in /( However we see in Figure 19 that there is no scaling with  cl , save for the slight drop-off at large cloud sizes due to the steeper scaling of the surface area there. Despite this, we find that  in still scales with cloud size, but this is because in our measurements,  in ∝  3/4 turb ∝  1/4 cl .Figure 20 shows this scaling, which is consistent with subsonic Kolmogorov turbulence.There are two possible interpretations here.The first is that our method of measuring  turb is sensitive to the scale of the cloud, and that the scaling comes from measuring velocity differences at the surface on cloud scale separations.The second, which we lean towards, is that that the turbulent length scale here is not set independently by each individual cloud, but is instead a common scale over all the clouds.This is motivated by the observation that all the clouds are evolving in the same common turbulent wind, where the turbulent velocity field has been driven on the outer scale of the  The remaining panels in Figure 19 hence assume that there is no explicit dependence of  in on  cl besides that from  turb .The middle panel shows the scaling of  in / 1/4 w , which is consistent with the expectation of scaling as  3/4 turb .Again, there is a drop-off at higher values of  turb , but this is similarly attributed to the higher surface area for the largest clouds that we divide by to get  in (leading to a slope shallower by a power of −1/4 since  turb scales with cloud size).The lowest panel of Figure 19 shows the scaling of  in / 3/4 turb with  w , where we see the 1/4 scaling with cooling time characteristic of turbulent mixing layer growth in the strongly cooling regime.
Finally, combining all of the above, we can put in some numbers and directly compare to the model presented above for the mass growth rate of clouds (equations ( 17), ( 18) and ( 21)).The orange dashed line in Figure 21 shows our model with the typical values in our wind:  3 ∼ 30,  ∼ 10 7 K,  turb ∼  ∼ 66 pc, and  crit,shear ∼ 5 pc.For  turb , we use the orange dashed line from Figure 20, so  turb ∼ 10( cl /10 pc) 1/3 km/s.Note there are no scaling factors here.We find that the model agrees extremely well with the estimated mass growth rates of the clouds in the simulation.Here, it is important to highlight some complications in the simulation data here that are not immediately obvious.Our model predicts an overall scaling of  ∝  but we find that there is a selection effect stemming from clouds in higher pressure environment being found in patches of the wind which are being compressed.Since these patches are smaller in scale, corresponding clouds are likewise smaller.This artificially makes the slope of  in Figure 21 visibly shallower at small cloud sizes.

Dependence of Turbulent Mixing on Wind Shear
Lastly, we look at the scaling between the turbulent velocity  turb in the mixing layer on the cloud surface with the relative shear velocity between the cloud and the surrounding wind  rel .As discussed above, plane parallel turbulent mixing layer simulations find that  turb ∝  rel , close to linear.However, cloud crushing simulations observe that turbulence persists even when the cloud is entrained.We find results consistent with both findings and well represented by equation ( 20).This is shown in Figure 22.At high shear (high  vel ), there is a strong scaling of  turb with  rel shown by the orange dashed line, consistent with  turb ∼  rel  rel with  rel ∼ 1/15.However, as  rel gets smaller,  turb levels off at roughly the sound speed of the cool gas, shown by the teal dashed line.This supports the picture that the KH instability is but one way of generating turbulence in the mixing layer, and that turbulence can continue to persist and drive mixing even when clouds become entrained in the wind (see discussion above equation ( 20)).

Non-thermal Pressure Support
Figure 23 illustrates the lack of thermal pressure balance between the clouds and the background wind and the contributions of the two main sources of non-thermal pressure support in the clouds: magnetic pressure and turbulent pressure.The top left panel shows a slice of thermal pressure normalized by the mean over the entire slice.The clouds appear generally as under-pressured regions relative to their surroundings.The top right panel shows this more clearly, where the orange dashed line shows the mean pressure of cool gas in the orange shaded region, and the red dashed line represents the mean pressure in the hot gas marked by the red shaded region.The strong dip in the intermediate region where the cooling time is short indicates that the mixing layer itself is not well resolved (Fielding et al. 2020).From left to right, the bottom row of histograms show ratios of (i) thermal pressures in the cloud and the surrounding wind, (ii) magnetic pressure and thermal pressure in the cloud, and (iii) turbulent pressure and thermal pressure in the cloud.From the upper right and lower left panels, we see that just looking at thermal pressures, clouds are under-pressured relative to their environment by a factor of 2. Magnetic pressure support is not large enough to make up the different, with the magnetic plasma beta  only being ∼ 4.This is despite having  ∼ 1 in the ISM prior to the onset of the SN driven wind.Instead, the missing pressure support is provided by turbulent pressure, which we define as  cl ⟨ 2 ⟩ cl .The turbulent pressure is roughly equal too the thermal pressure (which is equivalent to having a turbulent velocity equal to the sound speed of the cool gas).Hence, pressure support in the cloud is provided mostly in equal parts by thermal and turbulent pressure support, with a minor contribution from magnetic pressure.Taken together the total pressure of the clouds is, on average, equal to the hot wind pressure.

Cloud Origins: Passive Scalars
The passive scalars we have employed reveal some interesting points about the origins of these clouds.Figure 24 shows histograms of concentrations of passive scalars in our clouds which track the amount of cloud material that was originally (right before the first SN) cold ISM gas ( 0 ), cool ISM gas ( 1 ), or mass injected by SN events ( 4 ).These histograms are constructed from all the clouds in the catalog which spans a range of times as denoted in Figure 5.We do not see any time dependence of the histograms within this window.We find that most of the clouds are comprised of gas that was originally part of the cool ISM, supporting the fragmentation origin of the clouds, with only a small fraction originating from the cold ISM gas.It should also be noted that in both cases, the gas likely mixed and cooled to different temperatures, either in the wind or cloud, as evidenced by the other passive scalars ( 2 and  3 ) in Section 3. The close to unity values of  2 and  3 show that almost all gas had at different points in time been both cold and hot.On the left, we show a slice of thermal pressure normalized by the mean at 14 Myr.We see that that clouds appear generally under-pressured as compared to their surroundings.This is shown more clearly on the right, where the orange dashed line shows the mean thermal pressure of cool gas in the orange shaded region, and the red dashed line represents the mean thermal pressure in the hot gas marked by the red shaded region.Bottom Row: From left to right: histograms of the ratios of (i) the ratio of thermal pressures in the cloud and the surrounding wind, (ii) the ratio of magnetic pressure and thermal pressure in the cloud, and (iii) the ratio of turbulent pressure and thermal pressure in the cloud.These histograms are constructed from all the clouds in the catalog.

Cloud-Wind Alignment
Idealized setups involving clouds accelerated in a laminar wind or infalling under gravity often show an initial spherical cloud developing a pronounced head-tail structure, a morphological prediction that observations have widely been compared to (e.g., Putman et al. 2012).Clouds being accelerated in a highly turbulent wind exhibit far more complex structures, with the turbulence generating 'tails' that result in a rich diversity in cloud morphologies, as can be seen in the lower panel of Figure 15, where all clouds are pictured with the bulk flow in the upwards direction.Only a minority show similarities to classic head tail structures.In general, the head-tail paradigm serves only as a first order indicator of the direction of bulk velocities.We can quantify the alignment between clouds and a turbulent wind by looking at the angular separation between cloud and local wind velocities.This is shown in Figure 25, where we find that this angular separation follows an approximately exponential distribution with a scale angle of ∼ 25°.While most clouds are well-aligned with the wind, a handful are significantly misaligned.The existence of such misaligned clouds is consistent with observations of clouds in the Fermi Bubbles, but to our knowledge there have not been any detailed studies of such clouds (in contrast with aligned clouds, e.g., Di Teodoro et al. 2020;Noon et al. 2023).

Connection to Small Scale Simulations
It is exciting that the physical insights and analytic scalings derived from small scale idealized simulations translate well to larger scales.This was demonstrated in Tan et al. (2023) for individual infalling clouds, and now in this work for a population of clouds in turbulent magnetized outflows.In particular, these models build upon the work on the physics of interfacial turbulent radiative mixing layers (e.g., Ji et al. 2019b;Fielding et al. 2020;Tan et al. 2021), minimizing the need for empirical tuning of parameterized models.This success is promising for future use of physically informed subgrid models in larger scale simulations.Nonetheless, some nebulous points remain-we briefly discuss some of these issues, connections to other work, and limitations/caveats in this section.
What is the physical process that determines the initial cloud distribution?Inverse power law distributions, where the probability of a quantity taking on some value varies inversely with the power of that value, are widely observed in both the physical and social sciences, with a large existing body of empirical evidence supporting their existence (Newman 2005).They are commonly referred to as Zipf's law in the discrete case, or Pareto distributions in the continuous case.Their theoretical origins are, however, much less certain.In particular, it is difficult to explain their seeming universality.A ∼ −2 exponent, which indicates that there is a constant contribution from each logarithmic bin, is ubiquitous, even within astronomy (Guszejnov et al. 2018).Take, for example, the stellar initial mass function, which has a roughly  ()/ ∝  −2.35 scaling at higher masses (Krumholz 2014).This has been attributed to turbulent fragmentation (Padoan & Nordlund 2002;Hopkins 2013), where initially turbulent gas fragments into clumps.In this case, fragmentation is driven primarily by gravitational collapse.Gronke et al. (2022) find in simulations of turbulent multiphase boxes that the cool gas breaks up into droplets which follow a similar distribution.This setup is akin to a population of clouds fully entrained in a turbulent background wind.They argue that this power law, and in particular the −2 exponent, arises from competition between fragmentation and growth, where growth is provided from both a multiplicative source (cooling induced growth) as well as an additive source (coagulation of multiple smaller clumps with a larger clump).1 Fragmentation, on the other hand is driven solely by the turbulent field.Likewise, Fielding et al. ( 2023) find the same mass distribution of cool clumps in turbulent magnetized boxes where the multiphase medium is formed via thermal instability.While our setup differs from these more idealized ones, many features are similar.Instead of a single large cloud in a turbulent velocity field, we have an ISM that is broken up by an expanding hot superbubble.More work is needed to understand this process of turbulent driven fragmentation in detail.
In systems where the density contrast  is high (≥ 10 3 ), there is still remaining uncertainty around the cloud survival criterion.In this work, we have presented a criterion based on the idea that the characteristic cooling time of the mixing layers must be significantly shorter than the shear timescale that is motivated by and consistent with results of wind tunnel simulations at high  (Abruzzo et al. 2022a).This criterion predicts clouds must be larger than the Gronke et al. (2022) criteria by a factor of √  in order to grow, and is consistent with phenomenological criteria (e.g., Sparre et al. 2020;Li et al. 2020).However, there remains further investigation required to pin down the relevant physical processes that determine cloud survival in high  environments like multiphase outflows.We do not investigate cloud mass loss rates in this work.The methods we use to estimate cloud growth rates break down for clouds that are getting destroyed (see Appendix D).Incorporating cloud destruction rates into this model would improve the coupling between the phases, and is especially important in parameter regimes where clouds are not expected to grow.Similarly, this also requires more work on the small scale simulation side.
A final point on the topic of cloud survival, while not directly relevant to this work but important to keep in mind, is that Farber & Gronke (2022) show that the story is more complicated when  cloud < 10 4 K.This is mainly due to the shape of the cooling curve, which Abruzzo et al. (2022b) also showed plays an important role.However, the survival criteria they find is equivalent to the Gronke & Oh (2018) criterion when  cloud ≳ 10 4 K as it is in our simulations.
What about magnetic fields?This perennial question has been the subject of many wind tunnel simulations of cool clouds (e.g., Gregori et al. 1999;McCourt et al. 2015;Grønnow et al. 2018;Hidalgo-Pineda et al. 2023).Despite that, their effect remains unclear.For example, Gregori et al. (1999) find that destruction can be enhanced by more rapid acceleration while McCourt et al. (2015) find that they aid survival via magnetic draping inhibiting shear instabilities (see also Banda-Barragán et al. 2018;Ji et al. 2019b;Grønnow et al. 2022).Sparre et al. (2020) find some enhancement in cloud survival for  ∼ 10, while Li et al. (2020) find no effect for  ∼ 10 6 .Hidalgo-Pineda et al. ( 2023) find that for small  (∼ 1), much smaller clouds are able to survive.However,  is large in our winds and only significant in the disc-we hence conclude that it is unlikely that magnetic fields play a significant role in influencing cloud acceleration or survival in multiphase outflows.Magnetic fields don't seem to inhibit mass growth rates either (Gronke & Oh 2020;Hidalgo-Pineda et al. 2023), although morphologically clouds are reported to have very different filamentary structures as compared to their hydrodynamical counterparts (e.g., Tonnesen & Stone 2014;Gronke & Oh 2020;Jung et al. 2023).Our model does not include magnetic fields, supporting their lack of impact on growth rates.Understanding why this is so is an interesting avenue for future work.While we do not compare to a hydrodynamical run without magnetic fields, we do not observe clear filamentary structures-which we attribute to a combination of a weakly magnetized wind and turbulence.In general, turbulence in the wind is the main generator of complex morphologies seen in the clouds.

Connection to Galaxy/Cosmological Scale Simulations
Having demonstrated that much of the insight garnered from smallscale simulations translates to more realistic larger-scale systems we can now address how these processes fit into the overall landscape of galaxy formation and global-scale simulations.The impact of capturing these multi-scale multiphase effects has been seen in isolated galactic scale simulations with self-consistently generated multiphase winds, which find that properties of the hot wind including temperature, density, and pressure fall off slower than expected with distance, and also travel slower than single-phase adiabatic winds (Fielding et al. 2017b;Schneider et al. 2020).These effects are consistent with expectations of mass, momentum and energy exchange between cool clouds and the surrounding wind (Fielding & Bryan 2022).The total cool gas mass flux in Schneider et al. (2020) decreases with distance, suggesting that the cool gas is being destroyed and mass loading the hot phase, even for largest clouds.This may be understood by applying the  crit,shear criterion for cloud growth as opposed to  crit,cc .In addition, they find that cool clouds are underpressurized by up to a factor of 10.While we find that turbulent pressure support is significant, this cannot fully account for the large factor.Given that we also find that regions of phase space where the cooling time is shortest are the most underpressurized, this can likely be attributed to lower resolutions (Fielding et al. 2020).Additional detailed studies of multiphase winds at even higher resolutions (or with conditions in which it is easier to resolve  crit,shear ) will be helpful in shedding further light on how the multiphase interactions shape the overall evolution of the winds.
Recent large cosmological simulations also exhibit galactic winds with self-consistently generated cool phases.While most have insufficient resolution to resolve any clouds, certain zoom-in and high resolution simulations can do so, albeit marginally.For example, multiphase galactic winds are seen in TNG50 (Nelson et al. 2019) and FIRE-2 (Pandya et al. 2021).The latter characterized the multiphase nature of the outflow by analyzing the contribution to fluxes from each phase and found results broadly consistent with our findings and past tall-box ISM patch simulations (e.g., Fielding et al. 2018;Kim et al. 2020a).
Correctly capturing the multiphase nature of galactic winds is not only essentially for accurately modeling the winds themselves, but as recent work has shown, may also be essential for capturing the correct regulation of star formation and thus galaxy evolution.When winds are able to separate into multiple phases, the hot phase, which has very high specific energy, heats and stirs the CGM, which prevents new star forming material from entering the galaxy (Fielding et al. 2017a), and the cold phase ejects material directly out of the ISM.Standard galactic feedback models cannot capture the high specific energy phase in particular and as a result may be missing important regulation mechanisms, as was recently shown using regulator and semi-analytic models (Pandya et al. 2022;Carr et al. 2023), as well as isolated galaxy simulations (Smith et al. 2023).It is uncertain if cosmological simulations will ever achieve the resolutions necessary to properly capture the physical scales relevant to multiphase wind launching and dynamics.In which case, simulations and models such as those presented in this work are important in being able to bridge the gap towards achieving the capability to include realistic multiphase outflows using subgrid techniques.

Implications for Observations
A key missing piece in our understanding of galaxy evolution is the amount of mass and energy carried by winds from the ISM into the CGM and beyond.Most observations of galactic winds come from probes that are sensitive to cool gas with  ∼ 10 4 K (e.g., Heckman et al. 1990;Martin 1999;Rubin et al. 2014), although in some rare nearby cases the hot phase is observable in X-ray (e.g., Lopez et al. 2020Lopez et al. , 2023)).Translating from observed quantities to an inferred mass flux is a difficult and uncertain process, however in almost all cases the inferred mass flux is up to orders of magnitude smaller than what standard theoretical models predict.For example, McQuinn et al. (2019) analyzed a sample of 12 nearby starburst dwarfs from the STARBIRDS survey and found mass loading factors of 0.2−7 compared to the 1−100 predicted from simulations.Concas et al. (2022) studied ionized gas outflows in 141 main-sequence starforming galaxies at 1.2 <  < 2.6 from the KLEVER survey, and found that the ionized gas mass only accounts for less than 2 percent of predicted values.Resolving this profound conflict between theory and observations requires a better treatment and understanding of multiphase outflows on both sides.Here, we have shown that the nature of feedback is likely to be dramatically different from standard single-phase galactic feedback models, which in part helps to relieve this tension.We can, however, also use this multiphase wind picture to help refine our understanding and modeling of observations.
In order to make the most of galactic wind observations, particularly new and planned spatially resolved emission observations (e.g., Reichardt Chu et al. 2022), new observational modeling paradigms that take the multiphase nature of galactic winds into account are required.In particular, our finding is that the vast majority of the readily observable cool gas is in the form of clouds with a fairly wellunderstood size distribution and a relatively small volume-filling fraction.Furthermore, because we have shown the properties of these cool clouds are closely coupled to the energy containing hot phase, with the two phases shaping each other's properties, future multiphase models may be able to constrain not only the mass flux (cool phase) but also the energy flux (hot phase).
Our findings also provide insight into the nature of the CGM.Observations in the CGM suggest that the cool phase is not in pressure equilibrium with the hot phase (Werk et al. 2014).One way of accounting for this discrepancy is the addition of non-thermal components.As previously discussed, our results suggest that turbulence within clouds is a significant contributor to pressure support, while magnetic fields are a minor actor, although other candidates such as cosmic rays may also provide further cold phase pressure support (Butsky et al. 2020).

Further Considerations
Despite modelling and simulating a more realistic turbulent magnetized multiphase system, we ultimately still make simplifying assumptions.The following are some such limitations and caveats, many of which are each interesting enough in their own right to warrant exploration in future work.

Additional Physics
The most obvious and direct of these is the inclusion of additional underlying physics which were not incorporated into this work, but could possibly affect such systems by significantly modifying the dynamics or thermodynamics of the outflow.Detailed investigations and treatments of these processes are thus needed to accurately access their impact and importance, in order that they can be properly accounted for when modelling the behavior of multiphase winds.One such source of uncertainty that could potentially have a large impact is cosmic rays, which have been shown to be able to modify the properties of multiphase winds (Huang et al. 2022;Armillotta et al. 2022).They can provide an additional source of non-thermal pressure support and affect cooling (Butsky et al. 2020), or accelerate clouds directly (Wiener et al. 2019).Explicit viscosity (e.g., Li et al. 2020;Jennings & Li 2021) and thermal conduction (e.g., Brüggen & Scannapieco 2016;Li et al. 2020) have also not been included here (although we expect that in turbulent mixing layers, mixing is generally dominated by turbulent diffusion (Tan et al. 2021), explicit conduction can affect observables (Tan & Oh 2021)).The inclusion of non-equilibrium cooling as well as more sophisticated non-equilibrium chemical evolution models could potentially be important, especially for predictions of observables such as ion column densities (Sarkar et al. 2022).We also assume solar metallicity and abundances everywhere.This is clearly an oversimplification.When metallicities of clouds and their environment differ, the mixing can have significant implications for observables (Gritton et al. 2014;Heitsch et al. 2022).Dust survival or depletion in these systems is also important to understand, in particular at lower temperatures (Farber & Gronke 2022).Closely tied to this is self-shielding.While we assume that the whole box is optically thin, leading to clouds with  ∼ 10 4 K, more massive clouds with  HI > 10 17 cm −2 are likely able to self-shield and hence posses cold cores (Farber & Gronke 2022).In fact, molecular gas is observed both in the Milky Way (Di Teodoro et al. 2020;Noon et al. 2023) and other systems such as M82 (Walter et al. 2002).

Geometry
Besides being limited in terms of spatial extent, the nature of our setup lacks the correct geometry to track the outflowing wind further into the halo.While wind properties are often converged with resolution, local box simulations like this are limited by the Cartesian geometry of the setup which restrict the reliability of quantitative predictions of wind properties, which are often not converged with respect to box height instead.In particular, such a geometry which lacks the inverse squared distance scaling does not allow the adiabatic expansion and corresponding subsonic to supersonic transition of steady state winds that are a hallmark of predictions from analytic galactic wind models like Chevalier & Clegg (1985) (Martizzi et al. 2016;Fielding et al. 2017b).This in part motivated our focus on a cloud-centric analysis, rather than an extended look at associated wind properties such as mass and energy loading that are typically done in such setups.

Looking Forward
A natural next step is to apply and test cloud models further downstream in the wind.Analytic models of cloud growth/destruction as they are carried out further in the wind such as the framework outlined in Fielding & Bryan (2022) can easily be extended and applied to populations of clouds in the manner we have done here.While we have determined a physically motivated evolution and initial distribution, validating the time evolution of such a distribution coupled to an expanding wind is the natural next step (Anthony Chow et al., in prep).Additionally, the scope of the simulations can be expanded by exploring the effects of the different mechanisms listed above (such as cosmic rays) on the outflow.This will provide us with a more comprehensive understanding of the physical processes that drive multiphase galactic winds.Finally, understanding the impact of varying the properties of the galactic environment, such as gas surface density and metallicity, on the outflow will allow us to test the robustness of our results and determine the degree to which they are dependent on specific initial assumptions.In the long term, such models can inform subgrid approaches to modelling unresolved multiphase outflows in large scale galactic and cosmological simulations, allowing us to study the macro impacts on galaxy formation and evolution.

CONCLUSIONS
At the frontier of the field of galaxy formation and evolution lies the challenge to understand the multiphase nature of the environment and its implications.Galactic winds are a key component of these complex ecosystems and can be driven by stellar feedback channels.Observations reveal them to be common and multiphase in structure.Theoretical modeling and simulations of galactic winds have become increasingly sophisticated, aiming to reproduce and understand these outflows.In this work, we have built on turbulent radiative mixing layer theory and applied this to understanding and modelling the formation of multiphase outflows.We have run 3D MHD tall box patch simulations with a clustered SN driven wind outflow, with a focus on analyzing and modelling the properties of the resulting seeded cloud population.The main findings are summarized as follows: • Seeded by Fragmentation: During the breakout phase, the hot expanding bubble propagates through a multiphase ISM.The multiphase nature of the ISM tends to lead to asymmetric breakouts, and causes the outflows to fluctuate in power and direction.More importantly, it also leads to the fragmentation of the ISM during the breakout, which seeds the resulting hot outflowing wind with a population of cool clouds.We leave a detailed study of this process to future work.Consistent with this formation history, we find that clouds are mostly comprised of gas that was originally part of the cool ( ∼ 10 4 K) ISM, rather than from cold ( ≲ 10 2 K) or hot material ( ≳ 10 6 K).
• Turbulent Clouds & Winds: The uneven breakout and the presence of these clouds induces large scale turbulence in the wind.The clouds gradually get entrained via mixing induced accretion of momentum.The turbulent Mach number in the hot phase of the wind is ∼ 0.3, and magnetic pressure is extremely weak, meaning thermal pressure is dominant.However, clouds have turbulent Mach numbers ∼ 1 internally and at their surfaces where mixing occurs.
• Clouds Exhibit Complex Morphologies: This turbulent environment naturally leads to complex cloud morphologies that do not always conform to the head-tail description.In some cases, clouds can appear to be extremely misaligned with the bulk flow.
• Cloud Size Distribution: Cloud sizes are well described by a power law distribution of / ∝ −2.While this has been observed in previous works, we find that this scaling seems to hold even with the inclusion of magnetic fields in the ISM.The lower and upper scale cutoffs are consistent with estimates of the cloud survival radius and the scale height of the disc, respectively.We find that this scaling appears early on during the SN breakout stage-consistent with the source of clouds being the process of fragmentation of the ISM.
• Cloud Survival: The critical radius below which clouds can survive in a turbulent wind are consistent with  crit,shear which is given by equation ( 16) as proposed in Abruzzo et al. (2023).This criterion is motivated by a combination of turbulent box simulations from Gronke et al. (2022) and cloud simulations with high  in Abruzzo et al. (2022a), who proposed that the survival criteria is set by  cool,mix <  shear .
• Cloud Growth: By combining analytic models for the surface area to volume relationship of clouds ( cl ∝  5/6 ) and mixing layer theory for the mass inflow velocity  in , we can model the growth rate of clouds.We find that the model predictions are a good match to what we observe in the simulations, including predicted scalings and estimated .
• Role of Wind Shear: At high relative shear velocities between clouds and the surrounding wind, there is a strong relationship between the shear and the turbulent velocities which drive mixing and hence growth.However, as clouds get entrained, the turbulence becomes independent of the shear and is roughly the sound speed of the cool gas.The KHI is thus not the sole driver of turbulent mixing.Post entrainment, turbulence must instead be sustained by cooling or cooling-induced pulsations.
• Non-Thermal Pressure Support: Because turbulent velocities within these clouds are high (roughly the sound speed within the cloud), turbulent pressure provides as much support as thermal pressure.Surprisingly, even though the initial ISM had  ∼ 1, clouds are much more weakly magnetized, likely due to significant mixing with the hot high beta wind.These sources of non-thermal pressure support bring the clouds into total pressure equilibrium with their surroundings, despite having a factor of 2 lower thermal pressure.
In summary, we have shown that many of the physical insights and analytic scalings from idealized small scale simulations of radiative turbulent mixing layers and singular clouds in wind tunnels translate well to larger scale, more realistic turbulent magnetized winds.The multiphase component of these winds (a population of cool embedded clouds) can hence be well modelled, allowing for physics informed subgrid prescriptions which account for unresolved coupling between the various phases to be utilized in galactic and cosmological simulations where resolution limits are prohibitive.While refinements are required (e.g., the survival of molecular gas), moving forward, proper treatment of the small scales in this manner promises to pave the way towards tackling burning questions that remain regarding the role of multiphase feedback in galaxy formation and evolution.We initialize 3D boxes with the analytical profiles at  = 10 4 K and see how well hydrostatic equilibrium is maintained.The error is smaller at higher resolutions and is negligible.There is a density floor at 10 −6 cm −3 .

APPENDIX A: HYDROSTATIC EQUILIBRIUM TEST
Figure A1 shows a simple test setup in hydrostatic equilibrium for two boxes at different resolutions (low; below fiducial).The box is at  = 10 4 K throughout with a mid-plane density of 10 cm −3 .There is a density floor set at 10 −6 cm −3 .In both cases, the error over time gets smaller, and is consistently smaller in the higher resolution box.These tests provide simple sanity checks of the initial conditions.In practice, outflowing gas quickly becomes the dominant effect.has some shift), while the blue and red lines show early and late time values from the simulation.At late times, there is some decay in kinetic energy since we do not continuously drive.The turbulence also becomes uniform throughout the box as expected.

APPENDIX C: EXTENDED TOWNSEND ALGORITHM C1 Exact Integration Scheme
We first review the original Townsend algorithm (Townsend 2009).This is an algorithm for computing the change in temperature due to isochoric radiative cooling over some given time interval.It is based on an exact solution by noting that if one knows the cooling function, one can in principle simply integrate this function over the given time interval.The algorithm is as follows.We want to solve the following cooling equation: We can define the dimensionless Temporal Evolution Function (TEF) as follows: In principle, this only requires that 1/Λ() be analytically integrable.This is basically a normalized measure of the time taken to cool/heat from  ref to .We can then integrate the cooling equation over a timestep such that the integrated cooling function becomes: where  cool is the single point cooling time defined as and hence we can update the temperature as In practice, this is done in 3 steps -computing the TEF, evolving it over the required timestep, and then transforming back to the new updated temperature.In the next section, we compute the TEF for the two most useful cases, piecewise power laws and piecewise linear functions.

C2.1 Piecewise Power Laws
Physically realistic cooling function often come in the form of piecewise power laws that have been fitted to more complicated underlying models.The construction of the TEF and its inverse for piecewise power laws is given in the appendix of Townsend (2009), which are as follows.We assume the following functional form: ) for some temperature bin  k ≤  ≤  k+1 .The TEF is then The constraint that  () is continuous leads to the recurrence relation with  ref =  ( ref ) = 0.For cooling(heating),  ref is any temperature higher(lower) than the current temperature and we construct the TEF for decreasing(increasing) temperatures.Hence for heating, we express the above recurrence relation as  k+1 =  k + . . ., starting from  ref < .The inverse TEF is thus

C2.2 Piecewise Linear Function
We can likewise compute the TEF and its inverse for piecewise linear functions of the form where  is the slope in the temperature bin  k ≤  ≤  k+1 .The TEF is then: and and its inverse is:

C3 The Extended Algorithm
The main benefit of the Townsend Algorithm is naturally that it is based on an exact solution and hence not sensitive to errors associated with temporal resolution.This makes it the best choice when using a simple treatment of cooling/heating that is only a function of temperature and does not have a functional dependence on quantities such as ionic abundances.One caveat is that thermal and hydrodynamical evolution are decoupled over the course of a single timestep.
It is hence still recommended to include a further constraint on the timestep past the CFL condition that addresses this issue (for example requiring that the timestep be below some fraction of the cooling time).
In Townsend (2009), the algorithm was only outlined for cooling.Including heating is non-trivial because heating scales differently with density compared to cooling.This means that the net cooling/heating function becomes density dependent, which means that we can no longer optimize by pre-computing all TEFs.Another challenge is the behavior of power laws near equilibrium points.As far as we know, the only attempt to include heating and tackle this problem is in Zhu et al. (2017).In their implementation, they assume piecewise linear functions and only one equilibrium point.We generalize this further to generally work with piecewise power laws, as most cooling functions are represented as such, and an arbitrary number of equilibrium points.
When calculating the net cooling/heating at a given density, we assume that the resulting table itself represents a piecewise power law, instead of being a linear combination of two piecewise power laws.The exception is that for bins that have a zero crossing, we assume they are piecewise linear functions.One advantage of the piecewise functions above is that we can use linear functions only for bins with equilibrium temperatures to interpolate smoothly across the bin, where a power law formulation breaks down.This is fine as long as we dot our i's and cross our t's when computing the TEF and its inverse.Note that by construction, the TEF for any point with net zero cooling/heating will have an infinite TEF.Hence in practice one should set the cooling/heating to some non-zero tiny number with the appropriate sign.Furthermore, for cooling/heating, we use the next higher/lower bin as the reference temperature and only calculate the TEF down/up to the next equilibrium temperature, as we must always remain within a cooling/heating region if we begin there.
More specifically, we implement a new cooling module as a class that is initialized at the beginning of the simulation with some specified cooling and heating table.It is assumed that both cooling and heating are represented as piecewise power laws.The class gives the user the ability to call a cooling function for some input temperature, density and timestep and returns the new temperature in cgs units.The class also allows the user to query the single point cooling time at some given temperature and density, along with the minimum/maximum temperatures.For the cooling implementation, we use two preallocated scratch arrays, one to hold temperatures and one to hold the net cooling.At the start of the function, we first check that the temperatures are within the bounds of the provided table.If it is not, we return the lower/upper bound instead.We then populate the scratch arrays with a copy of the temperature table and the net cooling table for the given density.Next, we figure out which temperature bin we are in and check if there is net cooling or heating in that bin.If the net value is zero, we return the current temperature.We now case on whether there is net cooling or heating.For either case, the steps are similar, but the directions in which we compute quantities are opposite, along with the choice of some indexes.We first use the next bin over as the reference bin, but check for the edge case where we are in an equilibrium bin.In that case we modify the next bin to be at the equilibrium temperature and set the net cooling in that bin to be zero.If we assume in a bin that cooling and heating have the following power law forms: . (C14) We also flag the bin to be linear.We next calculate the TEF recursively downwards, until either we reach the bounds or we reach another equilibrium temperature, in which case we again modify the bin and flag it as a linear bin.We then compute the current TEF followed by the TEF after the timestep.Since the bins with equilibriums are linear, there is some finite value of  that if exceeded immediately returns the equilibrium temperature there.Lastly, we compute the inverse TEF, accounting for the fact that we might have crossed several bins.This gives us the new temperature.

C4 Testing
Figures C1 and C2 show a test of the above algorithm, where we have set up a simple cooling and heating curve with  = 2 and computed the calculate new temperature as a function of timestep for a range of starting temperatures.The black dashed lines in Fig. C2 mark the expected equilibrium temperatures.

APPENDIX D: ESTIMATING CLOUD GROWTH RATES
Given the large time interval between full data outputs, and the resulting difficulty in matching and tracking individual clouds across multiple snapshots, we would like to be able to estimate the mass growth rate of an individual cloud based solely on its properties from a single time snapshot alone.
Here we compare several methods of doing this, which we test by applying to a controlled setup containing only a single cloud and where the mass growth rate of the cloud is tracked explicitly with high time resolution.This provides a 'ground truth' which we can use to test our various estimators against.The setup consists of a cool cloud, initially at rest, falling under gravity in a constant hot background as detailed in Section 3 of Tan et al. (2023) (more specifically this setup corresponds to the Λ 0 = 30 and  cl = 300 pc run).At fiducial resolution here,  cl is resolved by ∼ 12 cells.We estimate the mass growth rate in three different ways.The first method is by estimating the total cooling luminosity  ≡ ∫ L.If we then assume that radiative cooling balances enthalpy flux (i.e., in the sub to transonic regimes; Ji et al. 2019b), we can estimate the mass growth rate as  Q ∼ /   hot .This also assumes that the bulk of the cooling luminosity comes from the mixed gas in the turbulent mixing layer and that the contribution from the cool and hot gas components is negligible.In Figure D1, we show how this stacks up against the are estimates-having multiple methods provides a cross check that allows us to make this estimate of  more reliably.
When measuring mass growth rates of the clouds in our main simulations, we thus use the mass flux through a temperature isosurface  surface , but normalized to match mass growth rates estimated using the cooling luminosity m Q .This choice allows for more consistency when discussing the scalings related to both the inflow velocity through this isosurface and its computed surface area.
The normalization is done by finding the constant of proportionality  iso discussed above.In other words, we calibrate a scaling factor  iso such that  ∼  iso  surface , where  iso ≡ ⟨  Q /  surface ⟩.

APPENDIX E: 'GROWTH' IN A FAILED WIND
It is interesting to present some results which show what happens when the superbubble breakout fails to generate a wind.In Figures E1  and E2, we repeat some of the analysis presented in the main text for one such simulation run so as to highlight several important differences.Figure E3 shows temperature slices at various times in this simulation run.After the superbubble breaks out briefly and vents the hot gas, the hot material quickly begins to cool and falls back to the disc instead of sustaining a hot outflow.Figure E4 shows various quantities such as outflow velocity, density, and pressure at 18 Myr.Figure E1 shows the same analysis of the cumulative size distribution as Figure 16.Here, the distribution appears slightly steeper for larger clouds and is more consistent with  (> ) ∝  −1.2 .We note that this simulation has a higher resolution of 1 pc.It is possible that this slightly steeper slope is related to the higher resolution, but more likely it is due to the brief and weak nature of the wind-these clouds show very different growth scalings which we discuss next.We find that while the surface area to volume scaling above still holds,  in in the cooling wind is independent of both cloud size and turbulent velocity, suggesting that the cool clouds are not growing via the same mechanism of turbulent mixing layers.In Figure E2, we show that instead,  in scales linearly with  w (i.e., inversely with the cooling time), consistent with the picture above that the growth is instead directly driven by large scale cooling and condensation of the wind (much like a cooling flow).This paper has been typeset from a T E X/L A T E X file prepared by the author.Figure E3.Slices of temperature at various times from a simulation where the breakout is not successful in launching a wind on the top side and is disrupted too early by cold gas in the ISM.This causes the breakout to cool and fall back towards the disc, a cycle which is repeated several times.

Figure 1 .
Figure 1.Volume rendering of the main simulation used in this work that demonstrates the presence of numerous cool clouds (shown in opaque blue) embedded within the hot outflow (shown in transparent orange) powered by the SNe (shown as the light source) exploding within the disc mid-plane.

Figure 3 .
Figure 3.The SNe rate as as function of the stellar population age.This consists of the core-collapse SNe rate (< 44 Myr) and the Ia rate (> 44 Myr).

Figure 4 .
Figure 4.The various stages of a superbubble breakout, driven by clustered SNe.As the first SNe go off, a hot bubble expands outwards, propagating through the multiphase ISM, leading to a complex fractal structure.As the superbubble eventually breaks out of the disc, the ISM fragments and breaks up, thus seeding the hot outflow with cool clouds.

Figure 5 .
Figure 5.The outwards mass flux over time at four different heights showing the initial breakout and outflow above the disc.This then reverses sides twice over the timespan of the SNe.Dashed line represents negative values, while the red shaded region highlights the time window after the SNe rate peaks in which we analyze wind embedded clouds in the next section.

Figure 6 .
Figure6.Slices of temperature at various times from our fiducial simulation where the breakout successfully launches a wind above the disc, despite the outflow being interrupted in the middle.

Figure 7 .
Figure7.Slices of temperature, outflow velocity, number density, pressure and magnetic plasma beta for a single time snapshot 14 Myr after the first SN.Clustered SNe drive a hot outflowing wind which contains a population of cool clouds in this wind.

Figure 8 .
Figure8.Slices of the velocities in the three axis-aligned directions.The vertical velocities  z reflect the bulk outflow velocity, while the velocities  x and  y parallel to the disc plane reflect the strong turbulence in the wind. 2

Figure 11 .Figure 12 .Figure 13 .
Figure 11.2D histogram of turbulent Mach number  turb at heights above 0.5 kpc estimated using a filtering method to account for bulk flows and temperature .The orange dashed line shows the average  turb (∼ 1) within the orange shaded region while the red dashed line shows the average  turb (∼ 0.25) in the red shaded region.These probe the cool clouds and the hot wind respectively.

Figure 14 .
Figure 14.Top panel: Histogram of gas cell temperatures from a single snapshot.Most of the gas is in the volume filling hot wind, with a peak at  = 10 4 K corresponding to cool clouds.The red dashed line represents our threshold value for identifying cloud material.Bottom Panel: From left to right: slices plots of temperature, cloud material, and identified clouds.

Figure 15 .
Figure 15.Bottom Panel: Meshes representing the temperature isosurface of a range of identified clouds.The cloud in the top left is shown in greater detail in the top panel.Top panel: Zoom-in of a a cloud surface colored by v  .The wind direction relative to the cloud is denoted by a red arrow.The inset shows the surface area weighted histogram of v  . turb is the standard deviation of this distribution.

Figure 16 .
Figure 16.Cumulative cloud volume distribution at various time slices.The pink region indicates where clouds are resolved by less than 16 cells.The grey region shows the expected power law scaling of −1.Grey dashed lines show expected lower (for survival) and upper bounds on cloud sizes.
,−21.4 ≡ Λ( mix )/(10 −21.4 erg cm 3 s −1 ), and  wind is the Mach number of the wind.Alternatively,Li et al. (2020) andSparre et al. (2020) lay out a survival criterion  cool,wind <  life , where  cool,wind is the cooling time of the hot background wind and  life = f  cc is the predicted cloud lifetime.f is an empirically calibrated function, given by f

Figure 17 .
Figure 17.Distributions of the Mach number of the wind relative to the cloud (top) and cloud overdensity relative to the wind (bottom).Dashed lines show population means.

Figure 18 .
Figure 18.Surface area of clouds plotted against their size.Most clouds follow a  5/6 or  5/2 cl scaling denoted by the dashed orange line, although this relation appears to steepen for a handful of the largest clouds.
For strong cooling regimes, Tan et al. (2021) calibrate the equations above to simulations of turbulent mixing layers and find that  in ≈ 10 km s −

Figure 19 .
Figure19.Scalings of the inflow velocity  in with various parameters-the cloud size, turbulent velocity, and pressure (as a proxy for the cooling time).Observed scalings appear to match theoretical expectations.

Figure 20 .
Figure 20.The turbulent velocity measured on the cloud surface scales with the cloud size in a manner that is consistent with subsonic Kolmogorov turbulence.

Figure 21 .Figure 22 .
Figure 21.Mass growth rates compared to the predicted scaling from the analytic model (orange dashed line), show good agreement.

Figure 23 .
Figure 23.Top Row:On the left, we show a slice of thermal pressure normalized by the mean at 14 Myr.We see that that clouds appear generally under-pressured as compared to their surroundings.This is shown more clearly on the right, where the orange dashed line shows the mean thermal pressure of cool gas in the orange shaded region, and the red dashed line represents the mean thermal pressure in the hot gas marked by the red shaded region.Bottom Row: From left to right: histograms of the ratios of (i) the ratio of thermal pressures in the cloud and the surrounding wind, (ii) the ratio of magnetic pressure and thermal pressure in the cloud, and (iii) the ratio of turbulent pressure and thermal pressure in the cloud.These histograms are constructed from all the clouds in the catalog.

Figure 24 .Figure 25 .
Figure 24.Histogram showing that most clouds (for the time window spanned, see Figure5) are comprised of gas that was originally part of the cool ISM pre-SN ( 1 ), with a small amount from the cold ISM ( 0 ).The contribution from SN injected mass ( 4 ) is negligible.

Figure A1 .
Figure A1.A test of the initial hydrostatic profile for two low resolution boxes.We initialize 3D boxes with the analytical profiles at  = 10 4 K and see how well hydrostatic equilibrium is maintained.The error is smaller at higher resolutions and is negligible.There is a density floor at 10 −6 cm −3 .

Figure B1 .
Figure B1.A test of constrained turbulence.Turbulence is injected at the start of the simulation into a uniform box with periodic boundary conditions.The orange dashed line marks the expected kinetic energy as a function of height in the box at the time of injection, while the blue and red lines show early and late time values from the simulation respectively.

Figure C1 .
Figure C1.The cooling and heating functions used in our test case, with  = 2.The grey points show the values in each bin.
temperature in a bin that transitions from heating to cooling or vice versa is eq =   Λ  Γ  1/(   −  )

Figure C2 .Figure D1 .
Figure C2.The computed heating/cooling over time for a range of starting temperatures.Black dashed lines show the expected equilibrium temperatures.
shows  Q /  surface as a function of cloud size for the clouds used in the analysis in Section 4, with  iso = 0.05 represented by the orange dashed line being the scaling factor we adopt.

Figure E1 .
Figure E1.Cumulative cloud volume distribution at various time slices from the failed wind run.The pink region indicates where clouds are resolved by fewer than 16 cells.Grey regions show the expected power law scaling of −1.Grey dashed lines show expected lower (for survival) and upper bounds on cloud sizes.

Figure E2 .
Figure E2.The inflow velocity here is shown to scale linearly with pressure in the simulation run which fails to generate a wind.