Strength in numbers: A multiphase wind model with multiple cloud populations

Galactic outflows have a multiphase nature making them challenging to model analytically. Many previous studies have tried to produce models that come closer to reality. In this work, we continue these efforts and describe the interaction of the hot wind fluid with multiple cold cloud populations, with their number density determined by different probability density functions. To do so, we introduced realistic cloud-wind interaction source terms and a time-varying cooling area. We find that the model reproduces well results from small-scale hydrodynamic simulations, but exhibits a general destructive behavior both for a single cloud population as well as multiple ones. We show that including multiple cloud populations can alter the evolution of the wind drastically. We also compare our model to observations and show that the differential acceleration of multiple clouds can lead to a non-negligible velocity `dispersion' relevant for down-the-barrel studies. Furthermore, we compute the emitted cooling surface brightness and find it generally too faint to explain observed Lyman-$\alpha$ halos.


INTRODUCTION
Galaxies are not static systems, but in contrast they are subject to complex physical processes, taking place inside them at different scales, that determine their evolution (Zhang 2018).A process of particular importance is the competition between gas flowing in and out of the galaxy.
The amount of gas and material inside the galaxy is determined by the exchange of gas between inflows and outflows, thus this process is essential for the fate of baryons in the galaxy (Bell et al. 2003), as well as for star formation and for regulating the black hole growth.Outflows have their own importance, mainly because they enrich the Circumgalactic and Intergalactic medium with their material (Tumlinson et al. 2017) and subsequently affect the star formation of their host galaxies as well as because they help us understand the mechanisms that powered them.
The main mechanism powering the outflows is thought to be feedback from star formation (Rupke 2018), or AGN (King & Pounds 2015).Using different probes, observations show that these outflows exhibit a vast range of temperatures.An example of different probes are neutral hydrogen (e.g Walter et al. (2002)) for  ∼ 10 − 10 2 K, photo-ionised metals, for  ∼ 10 4 − 10 5  (McKeith et al. 1995;Martin & Bouché 2009;Westmoquette et al. 2011) and X-rays for  ∼ 10 7  (e.g Strickland & Heckman (2009), for more details see Veilleux et al. (2005); Veilleux et al. (2020) and references therein).These observational findings lead us to the conclusion that the galactic winds are essentially multiphase.However the multiphase nature ★ E-mail: harrynikolis@gmail.com of the winds makes the theoretical modelling of this kind of system a challenging task.
That is because the cool gas is subject to processes that tend to destroy it such as hydrodynamical instabilites (Kelvin-Helmholtz, Rayleigh-Taylor Chandrasekhar (1961)).The destruction process of an initially static cold "cloud" with size  cloud destroyed by mixing with a hot galactic wind of velocity  wind is characterized by a timescale of  cc ∼  1/2  cloud / wind , with  the density contrast of the two phases, usually taking values of at least  ∼ 100.The characteristic time needed for the hot wind to accelerate the cloud, and become entrained by the wind is given by  drag ∼  cloud / wind , so the destruction time is shorter by a factor of  1/2 leading to the destruction of the cold phase (Klein et al. 1994).
A new approach to the problem is by including radiative cooling to the process of mixing, which leads to the formation of a radiative turbulent mixing layer.This mixing layer is an intermediate area between phases, with temperature at the geometric mean between the two components  mix ∼ √  ℎ   (Begelman & Fabian 1990).From the nature of the cooling curves, this intermediate temperature gas cools much faster than the cooling rate of each individual component.In this regime, a parameter space exists in which, mass can be condensed out of the hot phase, causing the cold gas to grow and eventually survive, with a lot of effort being done on characterizing the growth with different survival criteria (Marinacci et al. 2010;Armillotta et al. 2016;Gronke & Oh 2018;Gronke & Oh 2019;Li et al. 2020;Kanjilal et al. 2021;Sparre et al. 2020;Abruzzo et al. 2022b).
The presence of turbulence and cooling, together with the complexity of the hydrodynamic equations make the system too difficult to model theoretically.Most of our insights for phenomena taking place during the interactions of the different phases, come from an extensive study of small scale hydrodynamic simulations, with or without cooling (see for example Scannapieco & Brüggen 2015;Brüggen & Scannapieco 2016;Goldsmith & Pittard 2017;Schneider & Robertson 2017;Braspenning et al. 2023 andNaab &Ostriker (2017) and references therein).These simulations study the interactions for a vast set of idealized initial conditions and provide us intuition and dependencies that can be used in most general cases.The most quantitative results in terms of characterizing the mass transfer rates between the phases have come out from simulations focusing on the radiative turbulent mixing layers (Ji et al. 2019;Mandelker et al. 2020;Fielding et al. 2020;Tan et al. 2021).The results of these simulations have offered us a deeper understanding and precise scalings for the characteristic mixing velocity between the two phases in the RTMLs.Despite the focus on simulations for solving the "cloudcrushing problem", an analytical model describing the multiphase nature is still useful, as it will provide a straightforward way to test predictions and compare with observations.Various approaches to model theoretically the galactic outflows have been made in the past, at first by using only a hot fluid Chevalier & Clegg (1985), and then by using a radiative cooling component (Wang 1995;Silich et al. 2003Silich et al. , 2011;;Thompson et al. 2015;Bustard et al. 2016).Attempts have been made in the past to model the multiphase nature of these outflows (see e.g Cowie et al. (1981); Suchkov et al. (1996); Zhang et al. (2017).The most recent attempt has been by using a two component prescription, a hot fluid and a cold cloud interacting with each other, including radiative cooling and gravity by Fielding & Bryan (2022) henceforth referred to as FB22.
In FB22, the authors have developed the complete set of hydrodynamic equations governing the evolution both of the hot galactic wind, as well as the cold clouds embedded inside it, including terms that account for the interaction between the two components.Despite the number of equations involved in these highly complicated interactions, the evolution of the clouds is characterized by astounding simplicity, as the survival or destruction of the clouds is based solely on a  =  turb  cool,mix factor.The major finding of this model is that heavier clouds can survive and be entrained by the wind, while lighter get destroyed.This is in accordance with the most survival criteria, that connect cloud survival with the initial geometry of the cold clump.
We build upon that, and present a model which consists of a hot fluid, and an arbitrary number of cold clouds, interacting with the hot one.Each cold cloud is modelled as a clump population, with each population having a different initial mass, and the number density of each population is determined by a probability distribution.This paper is structured as follows: First we present the analytic model, with details given on the construction of the cloud-wind interaction source terms, an equation regarding the tail evolution, and the introduction of multiple cloud populations.Then we discuss about the results of the model, starting from various limiting cases.In the beginning, we solve the model for a homogeneous background, in order to compare it with hydrodynamic simulations.Then we proceed to a single cloud case, that enables us to understand the phenomena at play and compare it to recent models.After these two limiting cases, we solve the full model for different cloud populations and various distributions.Lastly, we calculate the the emitted luminosity and surface brightness along the line of sight, in order to compare the model with observations.The final part includes a discussion about implications of the model, different physical phenomena that came up and future directions.

MODEL
We will build the model step by step, with the goal of introducing a distribution of cloud populations interacting with the wind.
Starting from the initial conditions, we assume a Chevalier & Clegg (1985) model inside the galactic disc, which will provide the initial values for the wind quantities, when the clouds are introduced.
The clouds are introduced gradually just outside the galactic disc and each population interacts with the wind dynamically.Therefore, in order to describe the system, we need equations for the evolution of both the wind and the clouds.
For the wind perspective, we will use throughout the general steady-state fluid equations including gravity, heating and cooling: 1 Here, , ,  are the mass, momentum and energy source terms for the cloud wind interaction, and  =    and  cloud are the wind and cloud metallicity, respectively, and  esc is the escape velocity related to the gravitational potential, as introduced in FB22.
In our case these interaction terms are due to the presence of multiple cloud populations, and are obtained from results of smallscale hydrodynamic simulations, as will be seen in the next sections.

Inside the Galactic Disc
For the region  <  disc , with no gravity ( esc = 0) and heatingcooling (L = 0), assuming no clouds inside the galaxy and with uniform mass and energy injection ( ,  related only to the uniform injection and  = 0 because of the absence of clouds), the fluid equations reduce to: which is the well known Chevalier & Clegg (1985) model.Solving this model will provide the initial conditions for when the clouds are introduced, in the  >  disc region.
The initial conditions we use are the ones by FB22, thus: with SFR the star formation rate, and  mass,hot the hot mass loading factor  mass,hot =  hot /SFR.A similar term  mass,cold can be used when we introduce cloud populations.When we deal with a single cloud population, the  mass,cold is the cold mass loading factor of this population.However, in the presence of multiple cloud populations  mass,cold will describe the total cold mass flux of the system, which means that it will be equal to the sum of the loading factors of each population, multiplied with an appropriate statistical weight.Mainly, the statistical weight for each cold mass loading factor will be determined by the probability to find the different populations having a specific initial mass flux.
For the energy source term we can conclude for  <  disc with the same reasoning: where we used the energy loading factor  energy =   SN and assumed that all the energy is injected from supernovae that release energy  SN = 10 51 erg in the time an amount of a total of 100 ⊙ is formed.

Outside the Galactic Disc
Outside the galactic disc, we must use the full equations 1-3.The source terms now include the presence of the clouds and to find their form, we must focus on how the wind interacts with a single cloud exchanging mass, momentum and energy.

The cloud-wind interaction
From the conservation of mass, we can express the growth of the clouds, as mass flowing from the hot phase to the cold phase.Then, the rate of the mass flowing into the clouds can be written as (Gronke & Oh 2019): With the  cool the surface area of the cloud1 and  in the inflow velocity.From the formed RTML, a turbulence velocity  turb develops.
Using results from small-scale hydrodynamic simulations Tan et al. (2021) the turbulence  turb has the form: cool,cold is defined as the cooling  cool with temperature and metallicity fixed by the cloud values:  cool,cold =  cool ( =  cloud ,  =  cloud , ).For the cooling time, we use the cooling curve from Wiersma et al. (2009) as in FB22.We do not employ explicit heating but use a temperature floor instead mimicking the impact.
From the findings of Tan et al. (2021), the inflow velocity changes form, depending on whether the system is in the weak or strong cooling regime.The parameter controlling the two regimes is Damkohler number Da (Kuo & Acharya 2012), which is the ratio of the outer eddy turnover time to the cooling time Da =  turb / cool .The two regimes correspond to Da < 1 and Da > 1 for the weak and strong cooling respectively.Motivated by the scalings of Tan et al. (2021) At the turnover point between the two regimes, the inflow velocity is found to be  in ∼ 16 km s −1 (Tan et al. 2021).We choose a cloud temperature of  cloud = 2 × 10 4 K. 2  In order for the two regimes to smoothly connect at a cloud temperature  cloud = 2 × 10 4 K, the equations for  in are: for Da mix < 1 and for Da mix > 1: with M turb =  turb / s,c .The turbulent velocity  turb is saturated at high Mach numbers (Yang & Ji 2023), so when M =  rel  s,h > 1,  turb ∼ M 0 , thus for M > 1 we have: For the effective area of the cloud  cool , we assume that the length of the cloud expands in the direction of the flow and forms a tail, in accordance to hydro-dynamical simulations.This is where the cooling of the hot mass takes place, leading to cloud growth.We, thus, assume with  the length in the direction of the flow.In order to model the dynamics of the area of the cloud, we use the following qualitative arguments based on 'cloud crushing' simulations Marinacci et al. (2010); Armillotta et al. (2016); Armillotta et al. (2017); Gronke & Oh (2018): at early times, part of the cloud retains its spherical shape while the tail has a small width.Hence, in that phase the assumption of a cylindrical geometry is an overestimation of the true area.At later times, it is also important to recall that clouds that are destroyed do not have enough mass to supply the tail formation.Because of the above arguments we assume Here,  ′ describes the evolution of the tail.For this quantity we 2 Observations (McKeith et al. 1995;Westmoquette et al. 2009) trace the cold gas to temperatures ∼ 10 4 K.However, due to uncertainties in heating processes the exact choice of the cold gas temperature is arbitrary(see, e.g., Wiersma et al. 2009;Draine 2011).Because of the sharp dropoff in the cooling curve between  ∼ 10 4 K and 4 × 10 4 K, this range of temperatures is most likely and used in previous studies (e.g.Tan et al. 2021,FB22).One might argue that the choice of  = 2 × 10 4 K is somewhat higher than most of these studies who employ 10 4 K but firstly turbulent heating especially at the base of the winds can raise the cold gas temperatures above this, and more importantly, the exact choice does not play a crucial role for the evolution of the wind since the most relevant cooling time is the minimum cooling time ∼ 4 × 10 4 K (Tan et al. 2021;Farber & Gronke 2022).
assume that as the cloud moves in the hot medium, the tail is formed with a velocity equal to the relative velocity between the cloud and the medium.Thus we can write: We can translate this relation to a dependence to the galactocentric radius r, covered by the hot volume filling phase: This qualitative argument is a straightforward generalisation of the relation for the cooling area by FB22 in order to promote it to a dynamical variable.
The radius of the cloud can evolve according to a cylinder-like geometry (Huang et al. 2020): However, as we mentioned before at earlier times, the tail undergoes a rapid expansion with a small width, while part of the cloud retains its spherical shape.By assuming that the radius evolves according to Eq.21, we underestimate the radius of the cloud for the first few  cc and hence the cooling area.Because of the above, we only allow the radius to evolve adiabatically instead of shrinking: where we introduced  cloud,0 ≡  cloud ( = 0).
When the tail is not formed, there is no reason to deviate from spherical geometry while the clouds shrink, therefore With all the above arguments we have constructed the cloud growth term  grow , which is related to radiative cooling.The next step is to construct a term that is related to the cloud losing mass.This term will be related to the instabilities that tend to destroy the cloud.
For the mass loss term  loss , we use the results of the simulations from Scannapieco & Brüggen (2015), who find that clouds actually survive on longer timescales than expected from linear arguments (Klein et al. 1994): The fudge factor  is a free parameter in our model, which we set to  = 2 to match with hydrodynamic simulations.The total mass exchange rate is therefore With the above quantities specified, the evolution of the cloud quantities in terms of the volume covered by the expanding wind are derived by FB22 and are described by, starting with the mass of the cloud: The cloud velocity equation consists of three terms including the momentum transfer between the wind and the cloud, the drag-term and a gravity term.Thus, the velocity evolution is given by with  c the circular velocity related to gravity.The metallicity equation that govern the system, are derived using terms only related to the exchange of metallicity between the phases: Equations 26-28 govern the cloud evolution in the dynamical wind background.In order to complete the system of differential equations for the whole cloud-wind model, we have to see how these equations affect the wind macroscopically.

Wind related quantities
The next step is to zoom out from the 'microscopical' picture of the cloud-wind interaction, in order to connect the quantities derived in the last chapter with the source terms in the wind equations.We mainly focus on the final set of equations here; a more detailed derivation can be found in FB22.
The rate of how the wind gains or loses, in other words, the mass source term is related to how each individual cloud gives or drains mass from the wind, is given as the momentum source term can be derived with the exact same procedure, adding a drag term due to the relative motion of the cloud in the wind: with the drag-term being: and D drag is the drag coefficient.The energy term is similar.Defining: we can write: This is the form of the source terms in the wind equations, caused by the interaction with a single cloud population.We generalise these source terms, by adding a number of different cloud populations with different initial masses.We ignore cloud-cloud interactions in our model, so no cross-terms will appear in the equations bellow.Thus, we have for the total mass, momentum and energy: with the terms  cloud ,  cloud constructed in detail in FB22.Each cloud population "i" has a different initial mass, and its number density is affected by the choice of distribution.We use discretised versions of continuous probability distributions in the following way: with  ( 0,cloud,i ) the probability density function for a given distribution.We assume that the initial cold mass fraction is picked up from a probability distribution: with  cold,t the total cold mass fraction and  0,cloud,i the initial mass for a cloud of the population "".By defining the initial cloud number flux to be: The clouds are injected gradually between  disc and 1.33 disk with a powerlaw, which means: with  disc = 300 pc and  = 10.The choice of  can be tuned without significant changes.Thus, we have for the number density of the clouds: Now all the quantities in the equations are specified, so we can go on solving the model for different conditions.

Multiple cloud populations
In the general case of our model, the wind interacts simultaneously with different cloud populations, meaning groups of clouds that have initially different masses.The number of different cloud masses is arbitrary and is a free parameter in our model.The cold mass fraction for every cloud population comes along with a certain probability, and this probability is completely specified by the choice of the probability distribution the cloud masses follow.
We investigate three different examples: • a lognormal distribution: with  the mean of the distribution and  being free parameters of our model.
• a powerlaw distribution: with  = 2 and the range of the powerlaw distribution is a free parameter for the model (it is the range of the initial cloud masses we choose).Note that the proportionality factors do not need to be determined as the quantity is normalized by the sum of each value.
• a delta-distribution: with the singular cloud mass  0 being the only free parameter.This distribution is identical to the single cloud mass case and allows us to study the impact of a more realistic cloud mass distribution.

Homogeneous background
First, we will study the case with a constant background medium.While this is not realistic for a galactic wind, it will allow us to compare our analytic model with hydrodynamical 'cloud crushing' simulations (Gronke & Oh 2018;Li et al. 2020;Abruzzo et al. 2022b) The differential equations that govern this system are the evolving rates of the cloud mass  cloud , velocity  cloud , metallicity  cloud and the length of the cloud tail .So the equations involved are: 45) This limited case even though not so physically interesting, provides a test for how the model can reproduce the results from small scale hydrodynamical ('cloud crushing') simulations.In this scenario, the equations are influenced by how the tail evolves, especially in the early stages of the evolution.More specifically, without including a cut-off to the tail length (like Eq. ( 18)), the growth of the tail -even though the clouds initially lose mass -leads to a regrowth of the cloud.This is true even for clouds losing almost all their mass ( cloud < 0.05 cloud,init ) after several  cc .This behavior is moderated by imposing the minimum on the length evolution which accounts to the fact that the length cannot be growing when the cloud is losing mass.
It is worth mentioning that the use of the cut-off does not rule out the possibility of the cloud to start growing again even after losing a significant part of its mass.This is an effect that is commonly observed in simulations and is possible in our model especially in the multicloud case.The role of the cut-off has the physical meaning of regulating the evolution of the tail, when the mass is not enough to allow its formation, and is also essential in order to reproduce the results of simulations such as Gronke & Oh (2018).Furthermore, in order to have the best agreement with the runs of the above hydrodynamical simulations we impose a fudge factor  = 2 for the  loss term in Eq. ( 24).
From Fig. 1, we see that the solution of the equations match qualitatively well with the results of hydrodynamic simulations carried out by Gronke & Oh (2018).We see that the evolution behaviour between the model and the simulations is matched (regarding wether the cloud survives or not for different sets of initial conditions).Even though we do not have a one to one correspondence with the simulations, this is generally not the purpose of models like this.In this figure, we show the case of  = 100 and we present higher overdensities in Appendix B (cf. Fig. B1 and Fig. B2 for  = 300 and  = 1000, respectively).
Except for small disagreements (e.g., for  = 1000,  cool,mix / cc = 5), we find that our model reproduces the behaviour of the mass evolution (shown in the left panel) reasonably well.For the case of the velocity (right panel), we see that the values for the model are generally lower.This is expected, as we overestimate the mass evolution, and due to momentum exchange, this leads to an underestimation of  the velocity.We observe, for instance, that in our model, as well as in the simulations, clouds initially tend to lose mass and then start regrowing again, with a tendency to keep growing as the evolution continues.This behavior, as mentioned earlier changes continuously as we vary the parameter range, so it happens even in the cases that clouds lose almost all their mass rapidly.This effect is capped in the cloud-destruction case by the limit for the evolution of the tail, which forces the clouds that are instantly shredded to not be able to form a tail and eventually regrow.
Figure 2 shows the result of the model for a larger parameter range.Specifically, we show the mass after 30 cc as a function of  cool,mix / cc for various overdensities and Mach numbers.One can see that our model agrees with the survival criterion of Gronke & Oh (2018) within an order of magnitude.
There's a disagreement for larger overdensities where the models predicts mass growth even for  cool,mix / cc > 1, however, for  > 10 3 only very few hydrodynamical simulations have been run thus far (see, e.g., Abruzzo et al. 2022a) and how the survival ratio potentially changes for such overdensities is yet unclear.We, thus, fix  = 2 for the remaining parts of this study and generalize next to a dynamical wind background.

Single cloud population
We proceed to increase the complexity of the model by allowing cloud-wind interactions with only one cloud population.This is the immediate next step, before moving on to the full multi-cloud case.
The simplicity of the single cloud approach will allow us to understand better the behavior and the properties of the solutions.This will help us separate between the effects of individual cloud evolution and collective effects from the introduction of distributions in the next sections.Furthermore, we can compare the results of our model with recent work (FB22), as we will discuss in more details in Section 3.3.
While the single case is presented as a limiting case here, the results are actually included and used in the multicloud case of section 3.4, as they are equivalent with using a delta distribution.A systematic study of this distribution is motivated by the fact that the behavior of the individual bins from this section can guide our intuition to the more general solutions, and will be used as a comparison to the results of different distributions.
The model now includes the wind evolution, so the initial conditions for our equations are determined by the Chevalier & Clegg (1985) model inside the galactic disc, and eventually the loading factors  mass,hot ,  mass,cold ,  energy as in Eq. 8, as well as the value of the Star Formation Rate.The use of these initial conditions help us connect the model to the galaxy characteristics, making it easier to compare with observables.The clouds are introduced gradually outside the galactic disc.We move on to present the solution for the equations of the model in two different cases.
In Figure 3, we display a scenario with lower star formation rate and loading factors.This is a case where clouds are destroyed, and   there is mass loss in the whole integration interval.From the plots one can see that clouds are shredded instantly after they are introduced.
As it is expected from the cloud evolution, the lightest clouds are shredded faster.As we move on to higher initial masses, heavier clouds lose mass more slowly and end up with around half of their mass.Even the heaviest clouds in this case are not able to grow or even retain their mass.
The next case is presented in Fig. 4.Here we use larger SFR and  mass,hot .One would expect that the clouds would survive in this regime, because values for the initial conditions are higher (see FB22).In reality we still see mass loss, even though the model exhibits a different behavior.Lighter clouds initially gain mass rapidly, and then start to lose mass again, because the  loss ∝  rel dominates over the  grow ∝  in term.This change is related again to the scalings of the inflow velocity.For the heavier clouds, we observe a destructive behaviour, but there are not significant differences between the initial and final value of their mass .For the lightest cloud  cloud = 10 2  ⊙ , the integration stops because the wind has cooled down to the cloud temperature.This is a common feature for the high SFR cases, and only occurs for light clouds that accelerate rapidly.In general, we observe a significant difference in the behavior of the lightest clouds for the different initial conditions.Mainly there is the case exhibited before with rapid destruction and the one here, where clouds first grow and then lose mass.Heavier clouds do not follow this pattern, and their evolution is in general more steady.This leads us to check the model for the whole parameter range.In Figure 5, we check the final value for the mass of the clouds at a radius of  = 30kpc for different values of the survival criterion of Gronke & Oh (2018), and for a bigger parameter range.We see that we have a general cloud destruction behaviour.As explained heavier clouds usually end up with a mass near their initial mass, as their mass evolution is really slow.Light clouds are more sensitive to the choice of the initial conditions.The case that light clouds end the integration interval with mass higher than their initial value does not mean survival.As we have seen in individual cases, the loss term dominates at some point in the interval, or the wind fails to sustain them.From this plot, we also see that we are in general out of the allowed regions for the survival criterion.We can conclude, that at least for our model, there is no indication to assume that the survival criterion of a homogeneous background holds for a dynamical wind.
A conclusion from the above paragraphs is that cloud-wind interactions are present, we see that the behavior of the model changes.Now, the clouds have a tendency to eventually die in the integration interval for almost all the parameter range of the initial conditions.This is related to the saturation of the inflow velocity, which scales as ∼  s,h as we can see from Eq 15, while the loss term is always proportional to ∼  rel .

Comparison with Fielding & Bryan (2022)
In this single cloud case, we can compare our findings with similar models such as the one proposed by FB22.The models differ in the choice of the cloud source terms as well as the length evolution equation and the equation for the  cloud,loss term.
The changes in the cloud source terms eventually lead to a different inflow velocity  in which is the term affecting the mass growth equation.The form we propose in our model, and specifically the dependence of the velocity to different physical quantities (see Eq 13,14), is similar to the one from (FB22) .The difference between the two is by using  cool,cold instead of  cool,mix and by using the exact scalings from the hydro-simulations of Tan et al. (2021).The use of the same scalings leads to small differences in the prefactors.Both of these changes contribute slightly to the behavior of the inflow velocity.The main difference in the behavior of the function comes from the cut off and the saturation we introduced in Equation 16and Equation 15.This different behavior is what determines whether clouds get shattered or not between the two models.
The tail growth function has a more natural evolution, governed by a dynamical Equation 20.The tail in the compared model is assumed to be already formed when the clouds are introduced.This assumption gives higher values for the cooling area from the beginning, boosting the cloud growth term.
The terms related to loss for the cloud mass have similar functional form, but differ by pre-factors.The fiducial prefactors of FB22 might result in difference in at most one order of magnitude between the models, depending on the choice for the value of the prefactor.This effect individually contributes slightly to the general results, but is one of the factors for the difference between cloud survival or not between the models.Furthermore, our equation for loss includes the delayed destruction timescale of Scannapieco & Brüggen (2015).
Another difference that contributes to these results is that in the second model, mass evolution is determined by a  factor that differs by the  cool,mix  cc criterion for  > 100.This will lead to differences as the common initial conditions for the dynamical wind background usually gives values  > 1000.The mass evolution equations are built to depend directly on this  factor, which will lead to different survival criteria between the two models.In our case, there is no 'built in' factor that determines the evolution, with the  cool,mix / cc being satisfied rather naturally in the homogeneous case.
In Figure 6 we present the differences between the models for the cloud mass, inflow velocity and tail evolution for two different scenarios.The upper series of plots correspond to the instant cloud dying scenario, while the lower to a case with higher SFR.In the second case the two behaviours differ.We see that in a dying cloud scenario, the behaviour of the masses is similar.However, for a case with higher SFR,  M,hot , the second model achieves cloud survival.The length of the tail in our case evolves and asymptotically approaches similar values as the other model for higher radii.The inflow velocities in our model always start with higher values but quickly their values become lower.In our case, we can see the two different cappings for the inflow velocities as introduced in Equations 15,16.
It can now be seen clearly that the main differences between the models that lead to the different results are the different inflow velocities used.If we focus on the lower panels in Fig. 6, it is interesting to observe that mass loss in our model always starts near the point where the inflow velocity gets lower than the one from FB22.Because the most important difference between the two inflow velocities is the bounds used, the destructive behaviour in our case is connected to these terms.

Multiple cloud populations
After exploring the various limiting cases, we can now check how the model behaves in its full form, using multiple cloud populations, with the initial mass of each population determined by different probability distributions.
As mentioned in earlier sections, we use three different probability distributions, a log-normal, a power-law ( ∝ 1/ 2 cloud,init with a varying cut-off ranging from [ min , 10 7  ⊙ ].The values for the lower bound change from  min ∈ [10 2  ⊙ , 10 7  ⊙ ] depending on the mean mass of the distribution we choose.)and a delta distribution (which is the same as the single cloud case).The model can be run with an arbitrary number of discrete cloud populations  sampled from the respective distributions.We confirmed that increasing this number from our fiducial value of  ∼ 20 does not change the results significantly.Because in our model clouds do not interact with each other, the effect of the distribution affects the wind, as there is a sum of number densities for the clouds that is included in the equations, producing a collective effect on the wind quantities.This effect then back reacts on the clouds affecting their evolution.The effect of the distributions to the behaviour of the model depends on the initial conditions we choose.
Examining the total mass for all populations as a function of radius, leads to ambiguities as they move with different velocities in the wind medium.We check the total ratio of the cold mass flux in the end of the integration interval for different values of the mean mass initially.By doing this, we investigate the tendency of whether the total cold mass available in the system will grow or not, when we change the mean mass.In Figure 7, we explore two different parameter regimes, with the left plot corresponding to a case, where all clouds die.We observe that by increasing the mean initial mass of the cloud populations, the total cold mass flux tends to increase.Here a log-normal distribution seems to help the total cold mass lose mass more slowly, while the power-law and delta do not have significant differences.The similar behavior between these two is expected, as the leading contribution for the power law is from masses only near the mean value, exactly as in the -distribution.As we have also seen in the single cloud case, here different initial masses, have significantly different evolution regarding the way they lose mass.Even one order of magnitude difference in the initial mass leads to an accountable difference in the percentage of mass they lose finally.When we choose a lognormal distribution, we allow these differences to contribute to the final outcome, leading to the difference from the other two distributions.By making the lognormal wider ( = 1.5), the effect is even more enhanced, as expected as we allow more and more populations to contribute, and the mass loss becomes slower.In the right panel of Fig. 7, we show an intermediate regime where some clouds survive while others do not.We observe that for heavier mean mass initially, the flux ratio decreases.This can be explained by the fact that small clouds gain mass faster than the heavy ones.However, we see that we have again a slight increase for higher mean mass values.In this case almost all clouds gain more mass initially in the interval and then start to lose mass again.The general behaviour is that while the heaviest clouds gain mass with the lowest rate, they also lose it with the lowest rate.So in our case, when all clouds enter the losing phase, the heaviest clouds start to lose mass with the lowest rate, and end up the interval with more mass (normalised to their initial value) than the next lighter ones.The differences between the distributions here are not significant.This behavior is expected, by analysing the evolution of each population.Here the total differences between the final and initial mass are not big, with clouds ending up the integration with slightly more or slightly less of their initial mass.Because of this, even if we use distributions that allow a wide range of masses to contribute, the total effect will not be significantly different, because we do not have big variations in the evolution between the populations.
In Fig. 8, we pick a log-normal distribution with  = 1 and explore how the behavior changes for different initial conditions.The lowest curve corresponds to the log normal curve in the left panel of Fig. 7, while the highest to the one in the right panel.Here we can see exactly how with increasing the SFR and the  mass,hot we move from a case where all clouds are shredded from the beginning, to the case where at least the lighter clouds tend to grow initially and end the integration interval with higher masses even though they have the tendency to die.From this plot we also see the trends for the evolution of different masses.Lighter clouds have significant differences for different initial conditions from instantly shredded to initial growth, while the heavy ones have a very steady behavior which either tends to be the same or does not differ by orders of magnitude.This sensitivity of the lighter clouds to initial conditions is a clear demonstration of the argument that heavy clouds have a steady evolution in contrast to the lighter ones.
The results presented in these plots are not specific to the chosen initial conditions.Mainly, for a case with instant shredding, there is a big difference on the amount of mass each population retains.When moving to cases with initial growth, there are not big variations on the amount of mass the clouds gain or lose.
Starting from the second scenario, the collective effect produced is as having a single cloud with the same mean mass, so the choice of the distribution does not really affect the results.The reason for this behavior is that because each cloud mass does not have significant changes in its evolution, so by including more cloud masses in the interaction does not lead to a large effect on the wind and then back to the clouds.
On the other hand, we observe cases where the choice of distribution plays a significant role to the results of the model.This is the first scenario, where the evolution of each cloud population has a big difference between different initial cloud masses.So when we use distributions, we include contributions from different masses, thus resulting in an observed effect to the results.
In general, if we focus on the individual behaviour of one cloud of each population we again observe the heavier clouds have a steady evolution, while the lighter clouds have bigger changes in their evolution along the integration interval.What is different with the single cloud case is that sometimes, if a large range of initial cloud masses is included in the model, the rapid growth of small clouds is enhanced.Specifically, we observe cases where lighter clouds tend to regrow with the presence of heavier cloud populations.The contribution of this effect to the whole system depends on how far are the lightest populations from the mean value.Because of this, it is difficult to see the effects when total quantities are examined.Apart from this scenario, the lighter clouds either are shredded instantly, or gain mass rapidly.The loss term always dominates leading to destruction, even if they end up with higher mass than their initial at the end of the interval.So even if the behavior of each individual cloud might not be exactly the same as the single cloud case, we do not expect large differences when we examine mean values for the model, regarding the general behaviour.

Lyman-alpha Halo Surface Brightness
Due to the cooling processes taking place at the cloud-wind interface, we have Ly emissions from our model.We calculate surface brightness along the line of sight as well as the total luminosity emitted, to connect our model with observations.One can relate the mass transfer rate from the hot to the cold medium (of a given cloud population ) to a total luminosity via (e.g., Ji et al. 2019;Gronke et al. 2021 where   is the Boltzmann constant and  p is the mass of the proton with The total luminosity is therefore given by with being the emissivity.This implies that the surface brightness at an impact parameter  from the galactic center is Figure 9 shows that the luminosity of the halos has values ranging from ∼ 10 41 erg s −1 to ∼ 10 42 erg s −1 which agrees well with the observations (e.g., Steidel et al. 2010;Wisotzki et al. 2018).Luminosity increases with the increase of the star formation rate, which is expected as  grow also increases with SFR.We also observe that the luminosity decreases for higher values of  mass,hot .With the initial cold mass flux ratio fixed to  mass,cold = 0.15 this is expected because by increasing  mass,hot , we get to the regime where we have more hot mass in the beginning, comparing to the cold mass.By having less cold mass spread in the volume, there are less wind-cloud interactions happening, resulting to less cooling, so the luminosity drops.The whole contribution to the luminosity value comes from radii near the initial injection of the clouds.This can be seen by the fast drop on the surface brightness profiles as well, from which we can deduce the reasons behind it.In Fig. 10, we study the solution for the surface brightness profiles.The calculations done here are presented in physical units, because we want to examine the behaviour of the solution, together with the components involved in the formula.The surface brightness depends on the value of the growth term, and the total number density of the cloud populations.We observe that the values of  cloud ,  grow drop several orders of magnitude during the integration interval.From this behaviour, we can correlate the drop in the surface brightness profiles, as well as the contribution from the central bins to the luminosity, to the same drops in the  cloud ,  grow .
By switching to astronomical units and including a red-shift parameter to the surface brightness profiles, we can compare the results with observed Ly surface brightness profiles (Steidel et al. 2010;Wisotzki et al. 2018;Lujan Niemeyer et al. 2022).However because of the same drop in the surface brightness calculated in our model, the results do not match the observations which typically are much flatter; dropping around only two orders of magnitude in surface brightness levels in the same range of impact parameter (e.g best fit from (Wisotzki et al. 2018)).Thus, we observe a significant difference between the surface brightness profiles of the model and the observation fits.
In conclusion, while cooling, multiphase galactic winds can produce significant Ly emission, this is too centrally peaked in order to explain observed Ly halos.This strongly suggests that radiative transfer effects are at play leading to the observed Ly halo shapes (e.

DISCUSSION
In this section, we will analyze here some phenomena that arise with the introduction of multiple cloud populations.We will describe some characteristics and implications of the model, and discuss future directions.

Negative relative velocities
A new behaviour for the case of multiple populations, is that the heavy clouds influence the lightest clouds.As we can also see in Fig. 11, the lightest cloud populations in the distribution sometimes  get accelerated rapidly and they get velocities higher than the wind ones.
Even though all the equations are automatically zero when  rel = 0, the behaviour near zero is important in this case.For the cloud velocity to never exceed the wind value, the terms in the cloud velocity equation must get asymptotically to zero as  rel approaches zero.If this is not happening, then the light clouds can get accelerated beyond the wind velocity.In order to explain the phenomenon mathematically, one has to focus on the  grow appearing in the cloud velocity equation 27 due to momentum exchange.
Because of the backreaction of the clouds on the wind, the gradient of the hot gas velocity can be steeper than the dynamical timescale of a clump population, thus leading to an 'overshoot' of the velocity.This happens when  dyn,hot ∼ ℓ/Δ over a certain lengthscale ℓ is shorter than the acceleration time of the cloud -which is  drag ∼  cl / rel or  grow ≡  cloud /  grow for ram pressure or momentum transfer, respectively.
Physically, this effect can be easily understood: as the clouds carry more momentum, they have a finite deceleration time -akin to bullets flying through air.
In order to determine whether this effect is physical or not, it could be needed to analyze the implications of the pulsations for the system.Further work understanding the momentum exchange term of the two phases, might point to the correct direction.The total effect of the phenomenon to the whole system though is not big, because it appears for mass values far from the mean, where the number density of the population undergoing this acceleration is significantly lower.

Reasons for cloud destruction
A general feature of our model is that we observe clouds being destroyed in the whole range of the parameter regime.
The domination of the loss term ∝  rel over the growth term ∝  in , is happening mainly because of the saturation of the turbulent velocity  turb ∝  s,h , as explained in previous sections.The growth term in that sense, gets a rather complicated form with different cappings.On the other hand, the loss term follows the simple form from linear theory, with the extension from Brüggen & Scannapieco (2016).A point forward would be to investigate more the loss term through simulations, especially for high overdensities  ≳ 1000, in order to understand the nature of the loss term better.
In the same sense, a limitation of our work that might lead to the cloud destruction is the use of both  grow and  loss simultaneously.More precisely, both growth and loss equations have been tested in cases where growth or destruction dominates, respectively, but the presence of both of them together is an assumption.For example the growth equation is derived from an application of the Gauss law, We show the mean velocity (straight lines) and the root mean square velocity (dashed lines).The quantities are presented for two different distributions: a lognormal distribution (blue) and a powerlaw distribution (orange).We use SFR=15 ⊙ /yr,  mass,hot = 0.4,  mass,cold = 0.1.We choose a powerlaw ∝  −2 cloud .The lognormal distribution has an initial mass range of  cloud,init = 10 2 − 10 7  ⊙ and  = 1.0.Both distributions have initially same mean mass ∼ 10 5.5  ⊙ and one has to assume a one way flow of mass between the two components.Another approach for the model could be used, with each term dominating when there is survival or destruction only, but then the decision of which term dominates and when, is arbitrary.A model that could lead to these approximations physically is more ideal but far more challenging to construct.

Connection to observations
In section 3.5, we computed the surface brightness profile originating from cooling radiation because of the mixing and inflow of hot gas.The vast majority (≳ 90%) of this energy will be radiated away via Lyman- emission of neutral hydrogen.
As stated, we generally find that while the overall energy budget is non-negligible (∼ 10 42 erg s −1 ), the surface brightness dropoff is much steeper than observed in high- Ly halos (Steidel et al. 2011;Wisotzki et al. 2018).As we explore in section 3.5, this is due to the fact that both the temperature as well as the  grow drop (by several orders of magnitude; cf.Fig. 10).
For other emission lines, the fraction of cooling radiation is of course even lower.This implies that extended H emission observed in (local) galaxies are predominantly powered by recombinationas suggested by previous studies (e.g.McCarthy et al. 1987;Strickland et al. 2002;Hoopes et al. 2003).For instance, the H surface brightness levels of M82 are between 10 −14 erg s −1 cm −2 arcsec −2 and 10 −16 erg s −1 cm −2 arcsec −2 at distances up to ∼ 5 kpc from the nucleus (Yoshida et al. 2019).For this redshift ( ∼ 0) we find total surface brightness levels of 10 −14 erg s −1 cm −2 arcsec −2 to 10 −18 − 10 −19 erg s −1 cm −2 arcsec −2 , i.e, 2 and 3 orders of magnitude lower.Note that the H surface brightness will be approximately another two orders of magnitudes lower than the total computed here (e.g.Draine 2011;Bustard & Gronke 2022).Whether these recombination events are due to photoionization or shocks is beyond this study but discussed extensively in the literature (e.g.McCarthy et al. 1987).
Nevertheless, emission studies of galactic winds can help constrain our model further, e.g., by nailing down the clump size distribution.Here, we agnostically studied both lognormal as well as powerlaw (with −2 slope) distributions with the latter being supported by theory (Gronke et al. 2021;Tan & Fielding 2023) as well as commonly observed in the ISM (e.g.Salpeter 1955;Sadavoy et al. 2010).However, what the clump size distribution is in winds is still unclear but can be constrained with current and future observations of local winds.
Another way of observing outflows, is naturally through 'down the barrel' absorption (e.g., Martin 2005;Chisholm et al. 2015;Li et al. 2023).Here, the comparison of theory and observations crucially relies on the mapping between real and velocity space and thus on models such as the one presented here.One particular interesting fact of the multi-cloud model is the increased apparent 'velocity dispersion' due to the differential acceleration.That is, because differently sized clouds accelerate at different rates, at a fixed distance (and thus also as an integrated measure), there is not just a fixed (mean) velocity but also a potentially large dispersion around that mean.
We illustrate this in Fig. 12, where we show the mean as well as the 'root mean square' velocity of one particular model as a function of radius.Clearly, the variance of the velocities is non negligible -and crucially depends on the cloud distribution.We can see an order of magnitude difference between a log-normal (blue lines) and a powerlaw (orange lines) -in spite of the fact that their mean velocities (shown as solid lines in Fig. 12) are comparable.Interestingly this dispersion is rather large (∼ 100 km s −1 ) and, thus, potentially larger than usual turbulence in galactic winds (e.g., Schneider et al. 2020, found values of ≲ 30 km s −1 ).A similar value for turbulent winds (∼ 100 km s −1 ) was found by Tan & Fielding (2023).Of course this dispersion is not isotropic and thus not a turbulent component, however, this is impossible to disentangle from down-the-barrel observations alone.
In conclusion, our model can provide insights into observations by directly comparing (mock) observations to theoretical expectations.This is specifically useful for absorption line studies where our multi-cloud model can explain potentially large velocity dispersions required to explain observations (see, e.g., Li et al. 2023).

Caveats of the model
Semi-analytic models such as ours crucially rely on basic theory and simulation work to map out the viable parameter space and provide the scalings used.It is, therefore, also important to cross-check the model predictions with results of these previous studies -something we tried to achieve in section 3.1 by comparing a homogeneous background, single-cloud case of our model to 'cloud crushing' simulations.
However, this has not been done in the entire (realistic) parameter space for a lack of such basic studies.In particular, • most simulation work focus on the low  ∼ 100, low M regime (e.g., Murray et al. 1993;Gregori et al. 2000;Abruzzo et al. 2022b).However, the winds of our model usually show much larger overdensities ( ≳ 1000) and Mach numbers (throughout the evolution M ∼ 3 − 10).More simulation is needed there to calibrate semianalytical models in this regime.
• similarly, most cloud-wind simulations focus on a homogeneous background (exceptions include Gronke & Oh 2019 for outflowing and Heitsch & Putman 2009;Grønnow et al. 2017;Heitsch et al. 2021;Tan et al. 2023 for infalling clouds).However, clearly a realistic wind background is varying and more work is needed in this direction.
• while arguably the  in term is well understood from detailed simulations of turbulent mixing layers (Fielding et al. 2020;Tan et al. 2021), the (evolution of the) surface area of the cold gas is much less so.Our model includes a temporal variation of the surface area, however, this is clearly a simplification of the fragmentation and coagulation processes (Gronke & Oh 2022) occurring in multiphase systems.This is in particular true for clouds close to the destruction limit, i.e., where neither the growth nor the destruction term can be neglected -making essentially all of our clouds hard to model.
• a related point is the neglection of 'ensemble effects', that is, the influence of clouds (and cloud populations) directly upon each other.Our model includes the indirect (i.e., via the wind) influence of cloud populations but not the (arguably more important; Cowie et al. 1981;Alūzas et al. 2012;Poludnenko et al. 2002;Banda-Barragán et al. 2021;Tan & Fielding 2023) direct ones.This clearly requires more (simplified) simulation work to study these effects.
• other physical processes we neglected here are magnetic fields, cosmic rays, viscosity and thermal conduction.While, e.g., viscosity and thermal conduction might not affect directly  in (see Tan et al. 2021, for the inclusion of thermal conduction in turbulent mixing layer simulations), they can alter the surface area drastically and, thus,  grow (Brüggen et al. 2023;Brüggen & Scannapieco 2016;Scannapieco & Brüggen 2015).Similarly, magnetic fields can play a large role in cold gas survival (McCourt et al. 2015;Shin et al. 2008;Hidalgo-Pineda et al. 2023).
While our model provides a next step towards realism, and is useful to compare theory to observations, it is important to keep these caveats in mind.Overcoming them will require more theoretical and computational efforts of the community.

CONCLUSIONS
We present an analytical model for the interaction of the wind with multiple cloud populations extending previous work by FB22.In order to model the interactions between the phases as well as the time dependent area of the clouds, we used findings from small scale hydrodynamic simulations.Specifically, we implemented scalings found in high-resolution turbulent mixing layer simulations by Tan et al. (2021) and compared the final model to 'cloud crushing' simulations (Gronke & Oh 2018).Furthermore, we introduced probability distributions to describe the number density of each cloud population.Finally, we made an attempt to connect our model with observations, by calculating the emitted surface brightness profiles.
Our main findings are: • In the homogeneous case we can reproduce results from hydrodynamic simulations fairly well in the  ∼ 100 and with slight deviations in the high  ∼ 1000 case.For the whole parameter range of the initial conditions, the  cool,mix / cc survival criterion is satisfied except for high , M in the region between  cool,mix / cc ∈ [0, 10].
• For the single cloud population in an adiabatically expanding hot wind background, clouds have the tendency to die, because  loss ∝  rel dominates.For high SFR and  M,hot lighter clouds at first accelerate rapidly and gain mass but always start losing mass after some point in the evolution.Heavier clouds have a more steady evolution for the whole parameter range.
• The main difference compared to the FB22 model, when we use a single cloud population, is the saturation of the inflow velocity  in , together with a dynamical evolution for the length of the tail.
• For the multiple cloud case, the clouds have again a tendency to get destroyed.The choice of the distribution is important for certain values of the initial conditions, mainly the ones that we see significant differences in the evolution of each cloud population.For higher values of SFR and  mass,hot the differences are not that significant and the model is not very sensitive to the choice of cloud distribution.Curiously, in some cases the lightest clouds are "pushed" to negative relative velocities, mainly when there is a big range of initial cloud masses, that span several orders of magnitude.
• The cooling luminosity values are around 10 42 erg s −1 , while almost all of the contribution to this value comes from the central region.The reason for that is the drop in the  grow and  cloud term, after the first few kpc.The same effect can be seen in the surface brightness profiles.Due to this, there is a significant difference between emissions due to cooling and the observed Ly halos as well as measured H profiles indicating radiative transfer and photionization effects (cf.§ 3.5).
Our analytic multiphase galactic wind model captures the variety of cloud sizes a realistic galactic winds do possess.However, our work also shows that more work is needed in particular in constraining the source and sink terms via simplified, multiphase simulations.case is calculated at the injection point for the clouds ( ∼ 300pc).We proceed to check the behaviour of the criterion, if we calculate  cool,mix / cc , at larger radii, where the density of the clouds falls, and the conditions become more idealized.
In Figure A1, we calculate the cloud mass  cl,final ( = 30kpc)/ cl,init ( = 0.3kpc) and present it versus the  cool,mix / cc evaluated at three different radii: ( = 0.6kpc,  = 2kpc, and  = 5kpc).Because the model has a general destructive behavior, even the points with  cl,final / cl,init > 1 will eventually die.Thus, we would expect most values to be in the down and right quarter.By moving to larger values of the radius, the points are shifted to the right.We can see that the agreement becomes better for larger radii albeit with a large scatter.

APPENDIX B: DIFFERENT OVERDENSITIES
We compare the solutions of our model in the homogeneous case, with different simulations carried out by Gronke & Oh (2018), to show explicitly the disagreement that emerges for higher values of the overdensity . Figure B1, shows the cloud mass and relative velocity evolution for two different values of  cool,mix / cc = 0.5 and  cool,mix / cc = 5, for overdensity  = 300 and M = 1.5.We see that in this case, even though there are differences in the form of the cloud mass function, our model manages to reproduce the general behaviour of the simulations: Clouds grow for the  cool,mix / cc < 1 and die otherwise.
Figure B2 exhibits the case of overdensity  = 1000.Here four different values of the ratio are presented with different colours:  cool,mix / cc = 0.005, 0.05, 0.5, 5.In this high  case, we can see a slight disagreement of our model with simulations.More precisely, for the  cool,mix / cc = 5 case, our solution at first loses mass, but eventually regrows, while the simulated cloud loses mass and dies.Interestingly, the solution's behavior is similar to the  cool,mix / cc = 0.5 case.This disagreement is the same as the one discussed in Section 3.1, where in Figure 2, we see a discrepancy of our solutions and the criterion over the region  cool,mix / cc ∈ [1, 10].Apart from that case, we see that the rest of the solutions follow the behavior of the simulations and clouds grow, although the precise evolution of the clouds is different from Gronke & Oh (2018).
This paper has been typeset from a T E X/L A T E X file prepared by the author.

Figure 1 .
Figure 1.The time evolution of the cloud mass is presented in the left hand panel, while the time evolution of the relative velocity is shown in the right hand panel.The solutions from the model of the paper are presented in straight lines in comparison with the results of hydrodynamic simulations from Gronke & Oh (2018) which are shown in dashed.In order to use the same conditions as the hydrodynamic simulations, we use  = 100,  cloud = 4 • 10 4 K and wind Mach number M = 1.5.The Pressure of the wind is  wind /kb = 10 4 K cm −3 .Different initial values for  cool,mix / cc are presented with different colors.

Figure 2 .
Figure 2. A scatter plot presenting the final values of the cloud masses, normalised to the initial ones, as a function of the  cool,mix / cc criterion.The different colors correspond to three different  = 100, 300, 1000 and the symbols to 3 different wind Mach numbers M = 1.5, 2.25, 2.5.The allowed regions by the criterion are the upper left box and the lower right one.

Figure 3 .
Figure 3.A four figure panel for the evolution of basic quantities for different initial cloud masses in the single cloud case.The colors correspond to the initial cloud masses, which range from 10 2 (  ⊙ ) to 10 7 (  ⊙ ).The initial conditions are SFR = 7M ⊙ /yr,  mass,hot = 0.1, mass,cold = 0.1, energy = 1.

Figure 4 .
Figure 4.In this panel the initial masses range from 10 2 (  ⊙ ) to 10 7 (  ⊙ ).The initial conditions are SFR = 15M ⊙ /yr,  mass,hot = 0.4, mass,cold = 0.15, energy = 1.The integration ends sooner for the lightest cloud because the wind has cooled down to the cloud temperature.

Figure 5 .
Figure 5. Overview of different wind solutions showing the final cold mass versus  cool,mix / cc evaluated at the cloud injection point  = 300  (see Appendix A for different radii).The final cold mass is computed at  ∼ 30kpc.The initial cloud masses range from 10 0 .5( ⊙ ) to 10 7 (  ⊙ ) and every marker's size is proportional to the corresponding mass.The Star Formation Rate ranges from SFR=0.001(  ⊙ /yr) to 70(  ⊙ /yr) and the  mass,hot = [0.1,0.3, 0.5]

Figure 6 .
Figure 6.In this panel, we present the cloud mass, the length of the tail and the inflow velocity  in as a function of the radius.The straight lines correspond to our model and the dashed ones correspond to Fielding & Bryan (2022), while the different colors correspond to different initial cloud masses. 2 different cases are presented here, with the first row corresponding to SFR = 5M ⊙ /yr,  mass,hot = 0.1,  mass,cold = 0.1, and the second row to SFR = 10M ⊙ /yr,  mass,hot = 0.5,  mass,cold = 0.1.The integration ends sooner for the cloud with  cloud,init = 10 2  ⊙ in the second row, because the wind has cooled down to the cloud temperature.

Figure 7 .
Figure 7.We present the ratio of final cold mass flux over the initial value after the clouds are introduced as a function of the mean initial mass of the clouds populations.The different colours correspond to lognormal, powerlaw and delta distributions.The straight  = 1 and the dashed lines correspond to  = 0.75, 1.5 for the lognormal distribution.The left plot corresponds to the case where SFR = 5M ⊙ /yr,  mass,hot = 0.1,  mass,cold = 0.1, while the right has : SFR = 25M ⊙ /yr,  mass,hot = 0.5,  mass,cold = 0.15.

Figure 8 .
Figure 8.In this case, a lognormal distribution with  = 1 is picked.The behaviour of the flux ratio is again presented with the different colors corresponding to different star formation rates, while dashed lines correspond to different  mass,hot , while  mass,cold = 0.15 remains constant.The highest curve in this plot corresponds to the highest lognormal curve in Fig. 7.

Figure 9 .
Figure 9.The halo luminosity as a function of the Star formation rate, with different colors corresponding to different  mass,hot .The case studied here is a lognormal distribution with  = 1 and mean initial mass ∼ 10 5.5 M ⊙ /yr g., Chang et al. 2023; Byrohl & Nelson 2023) -maybe not unsurprising given the large Ly optical depth.

Figure 10 .
Figure 10.In this panel we present the Surface Brightness along the line of sight for different values of the impact factor , together with different quantities involved in the calculation of the Surface Brightness.The different colors correspond to different Star Formation Rates, while we set  mass,hot = 0.3,  mass,cold = 0.15.We use for the calculation a a lognormal distribution with  = 1 and mean initial mass ∼ 10 5.5 M ⊙ /yr.

Figure A1 .
Figure A1.In this panel, we present the model solutions for a cloud mass evaluated at  ∼ 30kpc and different values of  cool,mix / cc .In each plot,  cool,mix / cc is calculated at different radius:  = 0.6, 2, 5kpc respectively.The initial cloud masses range from 10 0 .5( ⊙ ) to 10 7 (  ⊙ ) and every marker's size is proportional to the corresponding mass.The Star Formation Rate ranges from SFR=0.001(  ⊙ /yr) to 70(  ⊙ /yr) and the  mass,hot = [0.1,0.3, 0.5]

Figure B1 .Figure B2 .
Figure B1.The time evolution of the cloud mass is presented in the left hand panel, while the time evolution of the relative velocity is presented in the right hand panel.The solutions from the model presented in this paper are presented in straight lines in comparison with the results of hydrodynamic simulations from Gronke & Oh (2018) which are presented in dashed.In order to use the same conditions as the hydrodynamic simulations, we use  = 300,  cloud = 4 • 10 4 K and wind Mach number M = 1.5.The Pressure of the wind is  wind /kb = 10 4 K/cm 3 .Different initial values for  cool,mix / cc are presented with different colors.