© The Authors 2026

1 Introduction

Galaxies eject material into their surroundings through powerful gas outflows, which can be powered by processes related either to star formation, such as supernovae (SNe) explosions and stellar winds, or to active galactic nuclei (AGNe). Multiphase outflows, from ionized to neutral and molecular gas phases have been extensively detected both in the local and in the distant Universe (see Veilleux et al. 2020, for a comprehensive review). From a theoretical point of view, a thermally driven hot wind that can accelerate cold material away from the galactic disk has been described by analytical and semi-analytical models with various levels of complexity, both for SN- and AGN-driven winds (e.g., Chevalier & Clegg 1985; King & Pounds 2015; Fielding & Bryan 2022; Thompson & Heckman 2024). The multiphase nature of outflows has also emerged from high-resolution hydrodynamical simulations of galactic winds (e.g., Gatto et al. 2017; Kim & Ostriker 2018; Kim et al. 2020b; Rathjen et al. 2021), where typically most of the energy is in the hot (T ≫ 10⁶ K) phase and most of the mass (at least at the wind base) is in the cold (T ≲ 10⁴ K) phase (see, e.g., Li & Bryan 2020), which may eventually return to the disk or evaporate in the hot wind (e.g., Schneider et al. 2020). However, due to the complexity of the processes at play, to the limitations of the observational data, and to the wide range of scales involved, detailed comparisons between theoretical models and observations are still missing.

A unique laboratory to observe and study multiphase galactic outflows is the Milky Way (MW) galactic center. The detection of the γ-ray emitting structures extending up to latitudes of about 55°, called Fermi Bubbles (Su et al. 2010), and more recently of their X-ray counterparts (eRosita Bubbles, Predehl et al. 2020), are a clear sign of an energetic outflow event from the central region of our Galaxy. Despite this evidence, whether such an event is due to the activity of the supermassive black hole (SMBH) at the center of the MW (e.g., Guo & Mathews 2012) or to central star formation activity (e.g., Crocker et al. 2015; Nguyen & Thompson 2022), which is particularly high (~0.1 M_⊙/yr) in a central ring with a radius of a few hundred parsecs known as the central molecular zone (CMZ; see Henshaw et al. 2023), is still a matter of debate (see review by Sarkar 2024, and references therein). Warm and cold gas exhibiting kinematics consistent with an outflow have also been detected in UV absorption (e.g., Fox et al. 2015; Bordoloi et al. 2017; Cashman et al. 2023), neutral hydrogen (HI) line emission at 21 cm (McClure-Griffiths et al. 2013; Di Teodoro et al. 2018), and down to molecular gas emission (Di Teodoro et al. 2020; Noon et al. 2023; Heyer et al. 2025) and absorption (Cashman et al. 2021).

The plethora of multiwavelength data of the Galactic center outflow gives us the opportunity to study in detail the interplay between the different gas phases of a wind. In this work, we focus on the cool phase seen in 21-cm emission close to the center of the MW (see Section 2), which accurately traces the gas kinematics and therefore can be used to test theoretical models of the Galactic wind dynamics. In recent years, these data have been interpreted using simple kinematic models of a cold outflow that has a constant velocity or accelerates with latitude (see Di Teodoro et al. 2018; Lockman et al. 2020). Our main goal is to follow up on such models, developing a more physically motivated scenario where a hot, fast wind powered by SN feedback entrains the cold clouds. While the physical properties and dynamics of the HI clouds close to the Galactic center have recently been studied using idealized high-resolution hydrodynamical simulations (Zhang & Li 2024; Zhang et al. 2025), in this paper we adopt semi-analytical models, with free parameters that we tune in order to reproduce the observed cloud kinematics. Although clearly simplified compared to simulations, this approach is more flexible and allows for comparison of observational data with model predictions that depend on key physical parameters.

The remainder of this paper is structured as follows. in Section 2 we summarize the HI surveys that we use in this work, and we report the main observational features that we aim to reproduce with our models. in Section 3, we describe the details of our multiphase wind models. In Section 4, we report the main findings of this work, which we obtained by comparing models and data through a Bayesian analysis. Finally, in Section 5, we discuss the comparison of our best-fit models with additional observational data and the main uncertainties of our framework, while in Section 6 we summarize our work, and we outline our conclusions.

Fig. 1 Population of HI high-velocity clouds from the surveys of McClure-Griffiths et al. (2013) and Di Teodoro et al. (2018), taken with the ATCA and the GBT. Left: galactic latitude-longitude map. The data-points show the positions of the clouds and are color-coded by their vLSR. The white mask shows the region surveyed by the ATCA and the GBT. Right: HI cloud latitudes as a function of their velocities. The kinematic pattern shows a signature of acceleration (see main text and Lockman et al. 2020). The points are color-coded by the cloud mass and the full 1D mass distribution is shown in the inset panel.

2 HI observations

The results of this work are based on two HI emission surveys of a region close to the Galactic center, presented in McClure-Griffiths et al. (2013) and Di Teodoro et al. (2018) and taken respectively with the Australia Telescope Compact Array (ATCA) and with the Green Bank Telescope (GBT).

Combined, the two surveys span a region within Galactic latitudes of |b| ≲ 10°, with the first reaching up to |b| ~ 5° and the second focused instead at 5° ≲ |b| ≲ 10°. The exact extent of the surveyed region is shown by the white mask in the left panel of Figure 1.

The main finding of these surveys is the discovery of a population of high-velocity HI clouds, with local standard of rest (LSR) velocities, v_LSR, up to ~400 kms⁻¹ and a kinematic pattern that is inconsistent with the Galactic rotation and likely associated with a Galactic nuclear outflow (see Di Teodoro et al. 2018). The total sample consists of 213 different components¹. in the left-hand side of Figure 1 we report the positions of the 213 clouds in Galactic longitude (l) and latitude, with the color scale showing the cloud v_LSR. Note that the absence of clouds within |b| < 2° is mainly due to the confusion of the velocity components with the emission of the HI Galactic disk. In addition, clouds with v_LSR < 50 kms⁻¹ for the ATCA sample and with v_LSR < 75 kms⁻¹ for the GBT sample in general cannot be easily distinguished from the HI emission from the local MW disk and are currently excluded from the catalogs. A more systematic search for clouds with low v_LSR will be presented in Yu et al. (in preparation).

The most evident property of the HI cloud kinematics is shown in the right-hand side of Figure 1: there is a clear sign of acceleration with increasing latitudes, as clouds closer to the disk tend to have lower LSR velocities. This feature was already identified by Lockman et al. (2020), who explained it with an outflow model that accelerates from a velocity of about 200 km s⁻¹ in the inner regions to a velocity of 330 km s⁻¹ at a distance of 2.5 kpc from the center. The aim of the current work is to expand on this finding by creating physical models of a multiphase wind where the cold clouds are accelerated by the fast-moving hot phase.

Additional properties of both samples, including cloud HI column densities, masses, and radii, are reported in McClure-Griffiths et al. (2013) and Di Teodoro et al. (2018). The clouds have typical column densities of up to a few times 10¹⁹ cm⁻² and radii that go from a few parsecs to about 50 pc. In the right panel of Figure 1 we have color coded the data points with the estimated cloud masses. While there does not seem to be a significant trend of the mass with cloud velocity or distance from the center, we note that the majority of the clouds with v_LSR ≳ 250 km s⁻¹ have relatively low masses ($m_{\mathrm{H}\mathrm{I},\mathrm{c}\mathrm{l}}\lesssim {10}^2{M}_{\odot }$). We also report in the inset in the same panel the one-dimensional mass distribution (the values have been updated with respect to those reported in Di Teodoro et al. 2018). The clouds have masses that range from a few solar masses up to 10⁵ M_⊙. We used the same mass distribution in our modeling (see Section 3.2).

3 Modeling the wind

In order to interpret the data reported in Section 2, we developed semi-analytical models of a multiphase wind of the MW, with the main assumption that the outflow is powered by SN feedback due to the intense star formation in the CMZ (but see also Section 5.2)². The idea is based on the steady-state wind model of Chevalier & Clegg (1985): these authors have shown that solving the hydrodynamic equations of the conservation of mass, momentum, and energy assuming constant injection rates of mass Ṁ_inj and energy Ė_inj within an inner radius r_inj leads to a fast-moving, high-temperature wind that quickly becomes supersonic outside of r_inj. The hot wind in turn can entrain and accelerate the HI cold clouds, potentially leading to the kinematic pattern observed through 21-cm data. We built our model on the more sophisticated framework developed by Fielding & Bryan (2022), hereafter FB22, which describes the steady-state structure of a multiphase galactic wind, with the cloud-wind interaction parametrized based on recent windtunnel simulations. The details of our modeling are described below.

3.1 Hot phase

With respect to the classical model of Chevalier & Clegg (1985), FB22 considered the effects of the galaxy gravitational potential Φ and of the radiative cooling and heating L on the hot SN-driven wind. They then assumed that the spherically symmetric wind subtends a solid angle Ω, which in our study we assume being equal to the full 4π. Assuming a smaller solid angle, and hence a biconical wind, would only change the hot gas density and pressure normalization by an order-unity factor and would therefore not affect our final findings (see also Section 5.2). Following FB22, the three conservation equations for the hot gas can be rewritten, taking into account the physical assumptions outlined above, to describe more explicitly the evolution of the hot wind velocity, v; density, p; and pressure, P, as a function of $$\begin{array}{rl}\frac{\mathrm{\partial }\log v}{\mathrm{\partial }\log r}& =\frac{1}{1-{\mathcal{M}}^{-2}}(\frac{2}{{\mathcal{M}}^2}-\frac{v_c^2}{v^2}-\frac{{\dot{\rho }}_{\mathrm{i}\mathrm{n}\mathrm{j}}/\rho }{v/r}\left(\frac{\gamma +1}{2}\right)\\ & -\gamma \frac{{\dot{p}}_{\mathrm{i}\mathrm{n}\mathrm{j}}}{\rho v}+(\gamma -1)\frac{{\dot{\mathcal{E}}}_{\mathrm{i}\mathrm{n}\mathrm{j}}-\mathcal{L}}{\rho v^2}+\frac{\gamma -1}{2}\frac{v_{\text{esc}}^2}{v^2}),\end{array}$$(1) $$\begin{array}{rl}\frac{\mathrm{\partial }\log \rho }{\mathrm{\partial }\log r}& =\frac{1}{1-{\mathcal{M}}^{-2}}(-2+\frac{v_c^2}{v^2}+\frac{{\dot{\rho }}_{\mathrm{i}\mathrm{n}\mathrm{j}}/\rho }{v/r}\left(\frac{\gamma +3}{2}\right)\\ & -\gamma \frac{{\dot{p}}_{\mathrm{i}\mathrm{n}\mathrm{j}}}{\rho v}+(\gamma -1)\frac{{\dot{\mathcal{E}}}_{\mathrm{i}\mathrm{n}\mathrm{j}}-\mathcal{L}}{\rho v^2}-\frac{1}{{\mathcal{M}}^2}+\frac{\gamma -1}{2}\frac{v_{\text{esc}}^2}{v^2}),\end{array}$$(2) $$\begin{array}{rl}\frac{\mathrm{\partial }\log P}{\mathrm{\partial }\log r}& =\frac{\gamma }{1-{\mathcal{M}}^{-2}}(-2+\frac{v_c^2}{v^2}\\ & +\frac{{\dot{\rho }}_{\mathrm{i}\mathrm{n}\mathrm{j}}/\rho }{v/r}(1+\frac{\gamma -1}{2}{\mathcal{M}}^2(1+\frac{v_{\text{esc}}^2}{v^2}))-\frac{{\dot{p}}_{\mathrm{i}\mathrm{n}\mathrm{j}}}{\rho v}\\ & +((\gamma -1){\mathcal{M}}^2)\frac{{\dot{\mathcal{E}}}_{\mathrm{i}\mathrm{n}\mathrm{j}}-\mathcal{L}-\dot{p}v}{\rho v^2}).\end{array}$$(3)

Here M = v/c_s is the Mach number, with $c_{\mathrm{s}}^2=\gamma P/\rho$ being the sound speed and γ = 5/3 the adiabatic index; $v_{\mathrm{c}}^2=rd\mathrm{\Phi }/dr$ is the circular velocity; and v_2sc = −2Φ is the escape velocity. In this model, within the injection radius r_inj, we have a constant mass density injection rate ${\dot{\rho }}_{\mathrm{i}\mathrm{n}\mathrm{j}}={\dot{M}}_{\mathrm{i}\mathrm{n}\mathrm{j}}/(4\pi r_{\mathrm{i}\mathrm{n}\mathrm{j}}^3/3)$ and an energy density injection rate in ${\dot{\mathcal{E}}}_{\mathrm{i}\mathrm{n}\mathrm{j}}={\dot{E}}_{\mathrm{i}\mathrm{n}\mathrm{j}}/(4\pi r_{\mathrm{i}\mathrm{n}\mathrm{j}}^3/3)$, with Ṁ_inj and Ė_inj being the mass and energy injection rate, respectively. For $r>r_{\mathrm{i}\mathrm{n}\mathrm{j}},{\dot{\rho }}_{\mathrm{i}\mathrm{n}\mathrm{j}}={\dot{\mathcal{E}}}_{\mathrm{i}\mathrm{n}\mathrm{j}}=0$. We also assumed no momentum injection, i.e., ṗ_inj = 0 everywhere (as in the original model from Chevalier & Clegg 1985).

Throughout this work, we assumed that the hot gas is subject to a logarithmic gravitational potential with a circular velocity v_c = 150 km s⁻¹. As we demonstrate later, this is a simplification that is slightly inconsistent with our treatment of the dynamics of the cold clouds, but the gravitational potential only marginally influences the properties of the hot wind (if the wind densities are not too high, which is the case for our best-fit models, see Figure 9) within 2 kpc from the disk (our region of interest), so our results do not depend on this choice. As for the cooling, we used the fiducial cooling tables used by FB22, who adopt the solar-metallicity cooling curves from Wiersma et al. (2009), which also include the effects of the extragalactic UV background at z = 0. We checked that a variation in the metallicity of the hot (and/or cold, see below) gas within a factor of 2 would not significantly affect the wind properties, and hence our final findings.

Finally, the injections of mass and energy within r_inj are related to the star formation rate, SFR, through the so-called mass and energy loading factors, η_M and η_E, such that Ṁ_inj = η_MSFR and Ė_inj = η_EE_SNSFR/(100 M_⊙), where E_SN = 10⁵¹ erg is the energy deposited by one supernova explosion and assuming that one supernova occurs every 100 M_⊙ formed. We assumed SFR = 0.1 M_⊙/yr, which is consistent with current observational estimates of the S FR in the CMZ (see Section 5.2). The injection radius and the two loading factors are the first three free parameters of our model (see Table 1, where we also summarize our choices for the fixed parameters described above).

In the left panels of Figure 2, we show the velocity (solid curves) and density (dashed curves) profiles of three different models where we assumed the same injection radius, r_inj = 100 pc, and three different choices for the mass and energy loading factors: η_E = 0.3, η_M = 0.1 (model A, top row); η_E = 0.3, η_M = 0.6 (model B, central row); and η_E = 0.5, η_M = 0.6 (model C, bottom row). One can see how in all cases the wind accelerates within r_inj to then reach a terminal supersonic velocity, while its density decreases substantially as it expands. As expected, increasing the mass loading factor implies that the wind becomes more massive and slower (the density increases by about an order of magnitude from model A to model B, while the terminal velocity goes from more than 1000 km s⁻¹ to about 600 km s⁻¹), while an increase in the energy loading factor increases the wind velocity, without significantly affecting the density profile (see the differences between model B and C). The hot gas profiles of Figure 2 are in line with the predictions of the more simplistic model of Chevalier & Clegg (1985), meaning that both the gravitational potential and the radiative losses have a minimal effect on the wind properties (especially for models A and C). In the next section we discuss how different hot gas models affect differently the properties and the dynamics of the cold HI clouds.

Fig. 2 Behavior of the multiphase wind for three different choices of models (see main text): model A (top row), model B (central row), and model C (bottom row). On the left, hot gas velocities (solid curves) and densities (dashed curves) are shown as a function of the distance, r. The right panels show instead the cold cloud orbits (solid curves), indicating in particular the effects of drag (dashed curves) and mixing (dash-dotted curves) on the original ballistic orbit (dotted curves). The symbols mark the positions of the modeled cloud at the beginning (circles), after 5 Myr (squares) and after 10 Myr (rhomboids).

3.2 Cold cloud dynamics

The main goal of this work is to reproduce the properties and especially the kinematics of the HI clouds presented in Section 2. For this purpose, we modeled the orbits of a population of clouds using a modified version of the python package GALPY (Bovy 2015), which includes not only the effects of the gravitational potential, but also of the interactions between the clouds and the hot gas phase described above, in a similar fashion to the works of Afruni et al. (2021, 2022) and following the framework of FB22. We assumed that the clouds have solar metallicity³; a temperature T_cl = 5000 K (see Di Teodoro et al. 2018); and that they are continuously pressure-confined by the hot gas, such that ρ_cl = μm_pP/k_BT_cl, with μ = 0.62. We note that this value of the mean molecular weight is appropriate only for the hot ionized gas, while for the cold phase it is μ_cold ≈ 1.2: we adopted this simplification following FB22, but a more accurate mean molecular weight should have a negligible impact on our findings.

3.2.1 Gravitational potential and multiphase interactions

We adopted one of the MW potentials available in galpy, the MWPotential2014. While we refer the reader to Bovy (2015) for additional details, in brief this potential has a circular velocity of 220 km s⁻¹ at 8 kpc from the center and it is composed of three components: (i) a spherical bulge that is modeled as a powerlaw density profile with an exponential cut-off; (ii) a disk that is described by a Miyamoto-Nagai potential (Miyamoto & Nagai 1975); and (iii) a dark matter halo that is described by a Navarro-Frenk-White (NFW, Navarro et al. 1996) profile. While this is not the most recent or accurate representation of the potential of our Galaxy, it is the most convenient for the purpose of this work as its simplicity significantly increases the computational speed of the final model with only marginal differences compared to a more sophisticated potential, especially for orbits in the inner Galaxy.

While the gravitational pull of the MW tends to decelerate the clouds ejected from the disk, the hot wind has the opposite effect of entraining and accelerating the clouds compared to a simple ballistic model (see Appendix A). We treated this entrainment by following the framework of FB22, where the acceleration of a single spherical cloud of mass m_cl is driven, once the gravitational force is taken into account, by two additional terms. The first one is due to the drag force exerted by the hot gas (see also e.g., Fraternali & Binney 2008; Marinacci et al. 2011): $${\dot{v}}_{\mathrm{c}\mathrm{l},\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}=\frac{1}{2}{C}_{\mathrm{d}\mathrm{r}\mathrm{a}\mathrm{g}}\rho (v-v_{\mathrm{c}\mathrm{l}}{)}^2\frac{A_{\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s}}}{m_{\mathrm{c}\mathrm{l}}},$$(4)

where $A_{\mathrm{c}\mathrm{r}\mathrm{o}\mathrm{s}\mathrm{s}}=\pi r_{\mathrm{c}\mathrm{l}}^2$ is the cloud cross section, with $r_{\mathrm{c}\mathrm{l}}=[3m_{\mathrm{c}\mathrm{l}}/(4\pi {\rho }_{\mathrm{c}\mathrm{l}}){]}^{1/3}$. We assume C_drag = 1/2 for consistency with FB22, who adopted this choice as it is applicable to flows with high Reynolds numbers.

The second term is instead due to hydrodynamical instabilities that cause mixing at the cloud interface between the low-velocity cold gas and the high-velocity hot gas: if the cooling time is sufficiently short, the mixed gas cools down and increases the mass of the cloud at a rate ṁ_cl,growth. This mass transfer from the hot wind into the cloud leads then to a transfer of momentum that significantly accelerates the cloud, so that (see FB22) $${\dot{v}}_{\mathrm{c}\mathrm{l},\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}=(v-v_{\mathrm{c}\mathrm{l}})\frac{{\dot{m}}_{\mathrm{c}\mathrm{l},\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}}{m_{\mathrm{c}\mathrm{l}}}.$$(5)

For an in-depth discussion on how the cloud growth rate is determined we refer to FB22 and references therein, while here we briefly summarize the basic concepts behind it, all based on results from recent numerical hydrodynamical investigations. As already mentioned above, the main idea is that at the interface (or mixing) layer between the cloud and the hot gas, the material can cool efficiently in a time τ_cool, equal to the cooling time of the mixed gas, which can be assumed (see Gronke & Oh 2018) to have an intermediate temperature $T_{\mathrm{m}\mathrm{i}\mathrm{x}}=\sqrt{T{T}_{\mathrm{c}\mathrm{l}}}$ and metal-licity (which in our case is always solar given that both phases have the same metallicity value). The growth rate is then given by (see, e.g., Fielding et al. 2020; Gronke & Oh 2020) $${\dot{m}}_{\mathrm{c}\mathrm{l},\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}=\rho A_{\mathrm{c}\mathrm{o}\mathrm{o}\mathrm{l}}{v}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}{\xi }^{\alpha }.$$(6)

Here, v_turb is the turbulent velocity in the mixing layer, which we approximated as 0.1(v - v_cl) following FB22 (the factor 0.1 comes from the simulations of Fielding et al. 2020 and Tan et al. 2021); A_cool is the cooling area of the cloud; and ξ = r_cl/(v_turbτ_cool) is the parameter that determines whether the system is in a slow-cooling (ξ < 1, for which α = 1/2) or a rapid-cooling regime (ξ > 1, for which α = 1/4). The cooling area does not coincide with the area of the spherical cloud, since the cloud tends to be elongated due to the interactions with the wind and to form a wake. Assuming that the evolved shape of the cloud is approximated by a cylinder with radius r_cl and height (v - v_cl)t_cc, where t_cc is the cloud crushing time, the cooling area is taken to be the lateral surface area of this cylinder, scaled by a correction factor, f_mix (corresponding to the parameter f_cool in FB22), which represents the efficiency of mixing layer formation (see also Section 5.3). By using the definition of t_cc = x^1/2r_cl/(v - v_cl), where $\chi ={\rho }_{\mathrm{c}\mathrm{l}}/\rho$ is the density contrast between the cold cloud and the hot medium, FB22 defined the cooling area as $$A_{\mathrm{c}\mathrm{o}\mathrm{o}\mathrm{l}}=f_{\mathrm{m}\mathrm{i}\mathrm{x}}4\pi r_{\mathrm{c}\mathrm{l}}^2{\chi }^{1/2},$$(7)

where f_mix is left free to vary as a free parameter of our model. Using the definitions outlined above, the final equation describing the cloud mass growth is given by $${\dot{m}}_{\mathrm{c}\mathrm{l},\mathrm{g}\mathrm{r}\mathrm{o}\mathrm{w}\mathrm{t}\mathrm{h}}=3f_{\mathrm{m}\mathrm{i}\mathrm{x}}\left(\frac{m_{\text{cl}}{v}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}}{{\chi }^{1/2}{r}_{\mathrm{c}\mathrm{l}}}\right){\xi }^{\alpha }.$$(8)

Finally, to express the overall mass evolution of the cloud, one needs to also take into account its mass losses. The same hydrodynamical instabilities (like Kelvin-Helmoltz instability) that lead to the creation of the mixing layers that might rapidly cool and increase the cloud mass and momentum, may also be responsible for shredding the cloud and making part (or all) of it evaporate into the wind. Following again FB22, the mass loss rate can be expressed as $${\dot{m}}_{\mathrm{c}\mathrm{l},\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}={\rho }_{\mathrm{c}\mathrm{l}}{f}_{\mathrm{m}\mathrm{i}\mathrm{x}}4\pi r_{\mathrm{c}\mathrm{l}}^2{v}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b},\mathrm{c}\mathrm{l}},$$(9)

where $f_{\mathrm{m}\mathrm{i}\mathrm{x}}4\pi r_{\mathrm{c}\mathrm{l}}^2$ is the area over which the cloud is mixed and v_turb,cl is the turbulent velocity responsible for mixing the cold cloud with the hot wind. Assuming that the turbulent energy density is constant between the cold and hot phases (Fielding et al. 2020), ${\rho }_{\mathrm{c}\mathrm{l}}{v}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b},\mathrm{c}\mathrm{l}}^2=\rho v_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}^2$, one can rewrite Equation (9) as $${\dot{m}}_{\text{cl,loss}}=3f_{\mathrm{m}\mathrm{i}\mathrm{x}}\left(\frac{m_{\text{cl}}{v}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}}{{\chi }^{1/2}{r}_{\text{cl}}}\right)$$(10)

so that the total cloud mass evolution is given by $${\dot{m}}_{\text{cl}}=3f_{\mathrm{m}\mathrm{i}\mathrm{x}}\left(\frac{m_{\text{cl}}{v}_{\mathrm{t}\mathrm{u}\mathrm{r}\mathrm{b}}}{{\chi }^{1/2}{r}_{\text{cl}}}\right)({\xi }^{\alpha }-1).$$(11)

In the right-hand panels of Figure 2, we show the effects of the forces described above on the cloud orbits, assuming the same three models analyzed in Section 3.1. For these examples, we assumed a cloud with a mass of 5 × 10² M_⊙, starting from a height z = 50 pc and a cylindrical radius R = r_inj = 100 pc with an initial vertical velocity of 100 km s⁻¹ and an initial tangential velocity equal to the circular velocity at the starting R. For models A and B, we adopted f_mix = 0.1 and for model C f_mix = 0.3. One can see how in all cases assuming the presence of the entraining hot wind tends to accelerate and push the clouds to much larger distances (solid curves) compared to a simple ballistic model (dotted curves). Moreover, the increase of momentum due to the mixing (Equation (5)) seems to be the dominant accelerating process and necessary to bring the clouds out to distances of at least 1 kpc from the disk (dash-dotted curves), while considering only the drag force (Equation (4)) would bring the clouds out to at most a few hundreds of parsecs (dashed curves). One can also appreciate the difference between the three models: in model C (bottom panel), the high hot gas densities and the higher f_mix lead the clouds to gain higher velocities and to reach slightly larger distances with respect to the other two models; in models A and B the clouds have instead similar final orbits, meaning that a denser, slower wind has an overall impact on the cloud evolution comparable to that of a faster, less dense one. In all the panels on the right-hand side of Figure 2, we also mark the position of the modeled cloud along its orbit at the beginning (before it is launched), at 5 Myr and at 10 Myr.

The cloud behavior can be better understood by looking at the evolution of the mass of the cloud in the different models, in Figure 3. In model A, even though the cloud is still accelerated as a consequence of momentum gain from the hot wind, it loses more mass to the hot medium than it gains from it at all radii, resulting in a net mass loss (ξ < 1 everywhere). In model B, the mass growth rate is, in the first part of the orbit, higher than the mass loss rate (ξ > 1), meaning that the cloud acquires mass from the hot wind. However, the latter is slower than the wind in model A, explaining the similar cloud dynamics between the two models. Finally, in model C the cold mass grows more efficiently than in the previous two models and the cloud acquires gas with velocities close to 1000 kms⁻¹, leading the cloud to a much higher acceleration.

We finally note that, throughout this work, we only considered the effects of the hot wind on the cloud orbits, but for simplicity we did not implement how the exchange of mass and momentum with the clouds affects in return the properties of the hot gas, which are instead fixed once the loading factors and the injection radius are chosen. This is an approximation of the original model of FB22 that we discuss more in detail in Section 5.3.

Fig. 3 Cloud mass evolution as a function of the distance from the center, for models A, B, and C (see Figure 2). In all cases, the cloud starts with a mass of 5 × 102 M⊙, but in model A it continuously loses its mass and slowly evaporates into the background medium due to hydrodynamical interactions with the hot wind, while in the other two models it initially increases its mass due to the cooling of the mixing layers.

3.2.2 Population of cold clouds

In our modeling, we assumed a constant flow of cold clouds that are ejected from the central region of the disk and whose motion outside of the MW disk is driven by the forces described in Section 3.2. To fully define an orbit, one needs to assume a starting position from which the cloud is ejected, the three components of the initial velocity, and a cloud mass. Given that our model is axisymmetric, we only need to choose an initial galactocentric radius, R_in, and height, z_in. We uniformly extracted an initial R_in between 20 pc and the injection radius, r_inj (see Section 3.1), and z_in between $\sqrt{r_{\mathrm{i}\mathrm{n}\mathrm{j}}^2-R_{\mathrm{i}\mathrm{n}}^2}$ and r_inj, so that the initial cloud population is effectively distributed on a thin shell outside r_inj and the clouds are not located inside the sonic region, where the hot wind is still strongly accelerating. We set the initial tangential component of the velocity, v_t,in, equal to the circular velocity at R_in, which we calculated using GALPY for the adopted MW potential. We then assumed that the initial kick velocity, v_kick, another free parameter of our model, is aligned with the vertical direction, so that v_z,in = v_kick and v_r,in = 0. We note that one could also assume a random direction for the kick velocity, but we verified a posteriori that this does not have any effect on the results presented in this work. The spherically symmetric hot wind tends indeed to strongly affect the cloud motion and the initial direction of v_kick has a minimal effect on the final cloud orbit. The kick velocity represents the initial velocity given to the cloud due to its interaction with the hot wind within the injection region, and it should be understood as the result of the cloud being initially swept up by the superbubble created by the SN explosions (which we do not capture in the model). Finally, the initial mass of the cloud is randomly extracted from the observed HI mass distribution, in order to have cloud masses consistent with the data. In particular, we truncated the observed distribution between 30 and 3000 M_⊙, to ensure that in the final mass distribution (given that the cloud mass can vary with time, see Equation (11)) there would not be clouds with masses lower or higher than those observed.

For a single model realization, we created 20 different orbits, and we followed the cloud motion for 10 Myr, with an additional stopping criterion in case the cloud gets to a distance r > 2.5 kpc, given that these distances are not covered by the observational data (Section 2). The orbit integration time was chosen to match the maximum cloud lifetime estimated by Lockman et al. (2020) with their accelerated outflow model. Indeed, given that we expected our physically motivated models to produce clouds with intrinsic velocities comparable to those predicted by the previous kinematic models (in order to be consistent with the observed cloud kinematics), we also expected comparable cloud lifetimes. We finally populated these orbits with N_cl clouds⁴ by placing them at different positions in the orbits assuming that they are evenly separated in time (see Afruni et al. 2021). Each cloud in the final distribution is uniquely described by (R_i, z_i,φ_i,v_r,i,v_z,i,v_t,i, m_i), where the azimuthal angle φ_i was drawn randomly between 0 and 2π. We performed this procedure above (positive z_i and v_z,i) and below (negative z_i and v_z,i) the plane of the disk.

3.3 Likelihood and priors

Once the cloud population was created, in order to compare with the observations we converted the cloud positions and velocities, initially in the reference system of the galaxy, to galactic longitudes and latitudes (l, b) and to line-of-sight velocities v_los (equivalent to v_LSR), using the in-built functions of GALPY⁵. We then took into account of the observational biases, by removing clouds that would fall outside of the GBT and the ATCA detectable positions and/or velocities and clouds that would overlap with the HI emission from the MW disk (see Figure 1 and Section 2).

We then defined our likelihood by directly comparing the final model and observed distributions of galactic latitudes, longitudes, and line-of-sight velocities, through $$\ln p_{\mathrm{t}\mathrm{o}\mathrm{t}}=\ln p_l+\ln p_b+\ln p_{v_{\mathrm{l}\mathrm{o}\mathrm{s}}},$$(12)

where p_x are the probability values of a Kolmogorov-Smirnov (KS) test performed between the observed and predicted cloud properties, i.e., longitude, latitude and line-of sight velocity, respectively. We used the likelihood expressed by Equation (12) to perform a Bayesian analysis across the 5-dimensional parameter space. As we have seen above, a model is indeed entirely defined by five free parameters: η_E and η_M define how much energy and mass is injected into the hot phase of the wind; r_inj defines both the radius within which the injection of mass and energy takes place and the initial location of the cold clouds; v_kick determines the initial ejection velocity of the cold clouds; and f_mix determines the strength of the mixing between the different phases of the outflow (see Table 1). We adopted a flat prior on the injection radius 0.05 < r_inj/kpc < 0.3, in order to consider radii that are roughly consistent with the extent of the CMZ (see Henshaw et al. 2023); a truncated Gaussian prior in log η_E between −0.8 and 0, centered in −0.4⁶; a flat prior in log η_M between −2 and 0.2; a flat prior in the kick velocity 0.3 < v_kick,₁₀₀ < 1.3 (where v_kick is expressed in units of 100 kms⁻¹); and a flat prior in the normalization of the mass exchange, −1.5 < log f_mix < 0.

Fig. 4 Corner plot showing the posterior distributions of the five free parameters of our models. The vertical lines mark the positions of the 2.5, 50 (solid line) and 97.5 percentiles of the one-dimensional posterior distributions.

4 Results

4.1 Bayesian analysis and comparison with the data

To perform the Bayesian fit, we adopted the nested sampling method (Skilling 2004, 2006), using the DYNESTY python package (Speagle 2020; Koposov et al. 2022). The posterior distributions of the five free parameters (see Section 3.3) are shown in Figure 4. While we discuss below the meaning of these bestfit values (Section 4.2), it is immediately evident how there is a well-defined region of the parameter space that better reproduces the data outlined in Section 2. The median values of the five posterior distributions are: r_inj = 0.23 kpc, η_E = 0.42, η_M = 0.14, v_kick = 59 km s⁻¹, f_mix = 0.32. The most prominent degeneracies are present between the hot gas mass loading factor and the kick velocity and the mixing normalization: less massive hot winds require cold clouds with larger initial v_kick and/or higher f_mix in order to produce outputs more consistent with the observations.

We first investigate the goodness of the best-fit models in reproducing the observed kinematics of the HI clouds. This can be appreciated in Figures 5 and 6, where we compare the data directly with the outputs of a model with the five free parameters fixed to the median values of the posterior distributions. In Figure 5, one can see the cloud line-of-sight velocity maps in galactic longitude and latitude for the best-fit model (top) and the data (bottom). Even though the distributions present some expected differences, we note that the general kinematic pattern is very well reproduced by our model, both in terms of velocity and of cloud location across the sky. The two distributions above the top panel additionally show the normalized distributions of longitudes for model (teal, dashed) and data (orange), whose comparison is used in the first term of the likelihood in Equation (12).

The consistency between model outputs and data is even more clear by the comparison in Figure 6 between the latitudes of the clouds as a function of their line-of-sight velocities (model in teal and data in orange). The modeled clouds show a sign of acceleration with increasing distance from the Galactic plane, very closely following the pattern detected in the HI data (see right panel of Figure 1 and Lockman et al. 2020). This is confirmed also by the comparison of the 1-d distributions of latitudes on the x-axis (corresponding to the second term of the likelihood, Equation (12)) and of line-of-sight velocities on the y-axis (third term in Equation (12)).

We can therefore conclude that we found a class of physically motivated models that lead to a cloud population in good agreement with the observations of Section 2, in terms of velocities, cloud locations and cloud masses (by construction, given that in the models we extract cloud masses from the observed mass distribution, see Section 3.2.2). In the next section, we explore more in-depth the physical properties of the multiphase wind predicted by our best-fit models.

Fig. 5 Comparison of the outputs of the best-fit models with the observational data. The cloud kinematics is shown as a function of Galactic longitude and latitude, with models on the top and data on the bottom. The one-dimensional distributions of the observed (orange, solid line) and the model (teal, dashed line) longitudes are reported on top.

Fig. 6 Comparison of the outputs of the best-fit models with the observational data. The model (teal) and data (orange) line-of-sight velocities are shown as a function of the Galactic latitude. The marginal onedimensional distributions are plotted along the two axes for both the data and the model.

4.2 Properties of the best-fit models

We first focus on the hot gas phase and specifically on the mass and energy loading factors of the hot wind. These are shown in red respectively in the top and bottom panels of Figure 7 as a function of the spherical galactocentric radius, r, and are given by $${\eta }_{\mathrm{M}}(r)=\frac{{\dot{M}}_{\mathrm{o}\mathrm{u}\mathrm{t}}(r)}{SFR}=\frac{4\pi r^2\rho (r)v(r)}{SFR},$$(13) $${\eta }_{\mathrm{E}}(r)=\frac{\dot{E}(r)}{E_{\mathrm{S}\mathrm{N}}SFR/(100M_{\odot })}=\frac{4\pi r^2\rho (r)v(r)[\frac{1}{2}{v}^2(r)+\frac{3}{2}{c}_{\mathrm{s}}^2(r)]}{E_{\mathrm{S}\mathrm{N}}SFR/(100M_{\odot })},$$(14)

where the two terms in Equation (14) take into account respectively the kinetic and thermal energy outflow rates.

The width of the bands in Figure 7 represents the 1 - σ interval of 200 models with parameters extracted from the posterior distributions shown in Figure 4. Both loading factors increase up to ~150-200 pc, corresponding to the injection radius (the median value is shown by the black dashed vertical line), within which we inject by construction all the energy and mass coming from the supernovae. At r_inj both factors therefore reach a maximum, with η_E being roughly equal to 0.4 and η_M ~ 0.1. Both loading factors remain then constant at larger radii: while this is expected by construction for the mass loading, the roughly constant η_E indicates that the radiative losses of the hot wind are negligible for this class of models. It is important to note that the values of the two loading factors are consistent with the findings of recent high-resolution hydrodynamical simulations of an ISM patch, like TIGRESS (e.g., Kim & Ostriker 2018; Kim et al. 2020b,a). In terms of star formation rate surface density, the Galactic center is most similar to the TIGRESS R2 model, with ∑_SFR ≈ 1 M_⊙ kpc⁻² yr⁻¹, for which they found an energy loading factor ηE ≈ 0.6 in the central regions and then decreasing with the distance from the disk. This justifies our choice of introducing a Gaussian prior at a value η_E ~ 0.4 in our Bayesian analysis. As for the hot gas mass loading factor, it has been found to be relatively independent of the central star formation rate and to be of the order of 10% (see Kim et al. 2020b; Schneider et al. 2020): we stress that our recovered value of η_M, which is only set by the comparison with the observations, is interestingly very similar to this value and in very good agreement with these theoretical findings.

We discuss in Section 5.1.1 the rest of the properties of the hot phase, but given the value of the loading factors, we can already conclude that the hot wind in our best-fit models is physically justified and consistent with expectations from state-of-the-art theoretical models.

In Figure 8, we report instead the intrinsic properties of the HI clouds predicted by our best-fit models. The different curves show the cloud properties along 20 different orbits that were created fixing the values of the five free parameters to the respective median values of the posterior distributions. The curves are colored based on the initial cloud mass, from cyan (lower mass clouds) to magenta (higher mass clouds). The actual orbits are shown in the upper left panel, where one can see that the clouds start from a thick shell close to the center of the MW (at about r_inj) and are then ejected and entrained by the hot wind at distances of 1.5-2 kpc from the center, consistent with the observations. Clouds with lower masses (≲10² M_⊙) tend to be more affected by the hot wind and are entrained up to larger distances with respect to more massive clouds (≳ 10³ M_⊙), which instead are generally not able to reach heights higher than about 0.8-1 kpc within 10 Myr. This effect can also be noted by looking at the cloud velocity components (bottom panels): while we can see that all the clouds get accelerated in the radial direction as they increase their distance from the center, the vertical velocities of high-mass clouds never exceed 100 km s⁻¹ and actually tend to decelerate from the initial kick velocity. We have seen in Section 2 (Figure 1) that in the data there is no strong dependence of the HI kinematics with the estimated cloud mass, but that most of the clouds with very high velocity have low masses, which is consistent with this emerging feature of our models. As for the tangential velocities (bottom right panel of Figure 8), the clouds are ejected by construction with the MW circular velocity at the respective initial galactocentric radius, R_in, going from about 70 to 90 km s⁻¹. Because of the gain of momentum from the hot wind (which is by construction non-rotating) the tangential velocities tend then to quickly drop and have values of at most 20 km s⁻¹ at r > 1 kpc. This feature is also consistent with the observed kinematic pattern of the HI cloud population, which does not show any sign of global rotation (see Section 2 and Di Teodoro et al. 2018).

The central top panel of Figure 8 shows the evolution of the cloud masses. One can see how, with this choice of parameters, all clouds tend to lose part of their mass into the hot wind, meaning that the mass loss term dominates in Equation (11) and that ξ < 1 at all radii. In particular, smaller clouds tend to be disrupted more rapidly and at 2 kpc have lost more than 70% of their initial mass. Therefore, although the clouds do accrete some high-velocity material, which contributes to their acceleration, they experience a net mass loss to the wind. This is an important result of our work: even though our framework allows the clouds to lose or gain mass, in order to reproduce the observed kinematics the clouds have to overall evaporate (i.e., lose mass due to stripping of material from the interactions with the background medium) into the hot wind as they move away from the disk. We therefore speculate that these clouds would likely not be able to reach distances from the MW disk larger than a few kiloparsecs.

We calculated the average mass lost by the clouds at a radius r, $f_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}(r)=⟨m_{\mathrm{c}\mathrm{l}}(r)⟩/⟨m_{\mathrm{c}\mathrm{l},\mathrm{i}\mathrm{n}}⟩$, in order to calculate the average cold gas mass loading factor as a function of radius, assuming that the N_cl clouds in our model have been ejected for 10 Myr at a constant rate: $${\eta }_{\mathrm{M},\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}(r)=\frac{N_{\mathrm{c}\mathrm{l}}⟨m_{\mathrm{c}\mathrm{l},\mathrm{i}\mathrm{n}}⟩f_{\mathrm{l}\mathrm{o}\mathrm{s}\mathrm{s}}(r)/(10\mathrm{M}\mathrm{y}\mathrm{r})}{SFR}.$$(15)

In a similar fashion, the energy loading factor of the cold gas was instead given by $${\eta }_{\mathrm{E},\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}(r)=\frac{{\eta }_{\mathrm{M},\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}(r)SFR[\frac{1}{2}{v}_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}^2(r)+\frac{3}{2}{c}_{\mathrm{s},\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}^2(r)]}{E_{\mathrm{S}\mathrm{N}}SFR/(100M_{\odot })},$$(16)

where $v_{\mathrm{c}\mathrm{o}\mathrm{l}\mathrm{d}}=\sqrt{v_R^2+v_z^2}$, which excludes the tangential motion that is not related to the wind. In order to get realistic loading factors, we chose a value of N_cl that would give us a final number of clouds, after taking into account of all the observational biases (see Section 3.3) consistent with the observed one (213). The loading factors are shown in blue in Figure 7, with the width of the bands representing the 1 - σ uncertainty calculated in the same way as for the hot gas. One can see that the cold gas mass loading factor is slightly larger than that of the hot phase, especially in the inner regions⁷, while it decreases at larger radii due to the cloud mass losses into the hot wind. The energy loading factor is instead almost two order of magnitudes lower than that of the hot wind. This picture is in agreement with high-resolution magnetohydrodynamical simulations of patches of ISM, where the cold outflow usually dominates in mass, with loading factors lower than unity (and decreasing with height) but generally higher than those of the hot phase, while the hot wind dominates the outflow energy budget (e.g., Li & Bryan 2020). However, we note that our calculation of the cold loading factors was based on the total number of ejected clouds over time and might not be directly comparable with the loading factors extracted from simulations (e.g., Kim et al. 2020b; Schneider et al. 2020), which come directly from the ratio between the mass flux (or outflow rate) and the star formation rate surface density (or star formation rate). Moreover, our value represents only a lower limit, given that it depends on the number of detected HI clouds: as already mentioned, the number, N_cl, used in Equation (15) was indeed chosen such that after taking into account the observational biases discussed above (see Section 3.3) the number of ‘detected’ clouds would be similar to the observed one (213); if many HI clouds are currently missed by the observations due to biases not taken into account in our analysis, then N_cl, and consequently the cold gas mass loading factor, would be higher.

Finally, the top-right panel of Figure 8 shows the radii of the clouds as a function of the distance, r, from the Galactic center. Even though the clouds are losing mass due to the interactions with the hot wind, they are increasing in size given that we assumed pressure equilibrium between the two gas phases and that the pressure of the hot gas drastically drops with increasing r. Therefore, r_cl goes from less than 10 pc close to the injection radius up to more than 100 pc for the most massive clouds at distances between 1 and 2 kpc. We note that these sizes are slightly larger than those estimated for the observed HI clouds, where the largest cloud radius is approximately 50 pc (see Di Teodoro et al. 2018). This could be explained by at least some of these clouds being out of pressure equilibrium. The distribution of timescales traveled by the clouds in the final cloud population of our best-fit models goes from 1 to 10 Myr⁸, with a peak around 5 Myr (consistent with the previous estimate from Lockman et al. 2020), while the cloud crossing time is between 2-15 Myr, meaning that some of the clouds could still not have reached an equilibrium state. in addition to this, the discrepancy between the observed and modeled cloud sizes might also be due to uncertainties in the density profile of the wind, which we discuss more in detail in Sections 5.2 and 5.3.

Fig. 7 Mass (top) and energy (bottom) loading factors as a function of the galactocentric radius, r, for a class of best-fit models (see Section 4.2 for more details), for the hot gas (red) and the cold clouds (blue). The vertical dashed line marks the location of the median injection radius.

Fig. 8 Intrinsic properties of the HI clouds predicted by our best-fit models. Top, from left to right: cloud orbits, cloud masses as a function of r, and cloud radii as a function of r. Bottom: three components of the cloud velocities. The different curves are colored based on the initial cloud mass, from lower masses in cyan to higher masses in magenta.

5 Discussion

We have found in Section 4 that there is a class of physically motivated models, consistent with various theoretical expectations, that predict HI cloud properties and kinematics in good agreement with the available constraints from 21-cm data of a population of high-velocity clouds close to the Galactic center. in this section, we explore the comparison of our model predictions with additional observational data of the multiphase MW wind (Section 5.1) and we discuss more in detail the assumptions and uncertainties of our modeling (Sections 5.2 and 5.3).

5.1 Comparison with multiwavelength observations

5.1.1 X-ray data

Although the main aim of this work is to infer and interpret the dynamics of the cold HI clouds, the hot phase of the wind is a fundamental ingredient of our physical picture. The motion of the clouds is strongly influenced by the interaction with the hot gas and, without accounting for this phase, it would not be possible to accurately reproduce the observations (see Appendix A for more details). We have already discussed above (Section 4.2) that the hot diffuse wind has mass and energy loading factors that are in agreement with theoretical expectations. Here, we compare its properties with the observational constraints currently available from X-ray data.

In the three panels of Figure 9 we show the intrinsic properties (velocity, density, and temperature) of the best-fit hot wind predicted by our Bayesian analysis, as a function of the radius, r. The pink curves show the outputs of 200 models with parameters taken from the posterior distributions shown in Figure 4 between the 16^th and 84^th percentiles, while the red bands show the 1 - σ scatter at each radius. One can see that within the injection radius the wind is still accelerating and has densities of about 5 × 10⁻³ cm⁻³ and temperatures between 10⁷ and 10⁸ K. At the injection radius, the wind reaches then a terminal velocity of about 1500 kms⁻¹ and its density and temperature start to drop to reach values of about 10⁻⁵ cm⁻³ and 10⁶ K at r ≈ 2 kpc.

While current X-ray observations do not have the necessary spectral resolution to infer hot gas velocities, we can already compare our densities and temperatures with observational estimates of the hot material detected in the vicinity of the galactic center. A few years ago, using XMM-Newton, Ponti et al. (2019) discovered X-ray structures extended up to hundred-parsec scales from the center, which they named Galactic Center Chimneys (see also Nakashima et al. 2019). These are almost cylindrical hot plasma structures that extend more than 1 degree in latitude from the Galactic disk; they could be the channel of energy transport from the Galactic center to the Fermi/eRosita bubbles and might be powered by either SN explosions or activity of Sagittarius A^*.

Given their properties, these chimneys could potentially be related with the inner part of the hot wind predicted by our model. In the second and third panels of Figure 9, we report in black the observational estimates of the chimneys’ temperature and densities, which span a range in galactocentric radii from about 70 to 200 pc (see Ponti et al. 2019). Our models predict a hot gas phase that is on average hotter and less dense than these observations⁹. This might imply that with our models we are constraining a different component, which is still too faint to be detected by current X-ray facilities.

To further investigate the discrepancy between our models and the data discussed above, we ran an additional Bayesian fit where we fixed the mass loading factor to values close to 1 (imposing an additional Gaussian prior in η_M), in order to force the wind to higher densities, more compatible with the data from Ponti et al. (2019). The median values of the resulting posterior distributions are: r_inj = 0.23 kpc, η_E = 0.54, η_M = 0.98, v_kick = 44 km s⁻¹, f_mix = 0.22 (see also Figure B.1). The hot gas outputs of such model are shown by the dashed curves in Figure 9: the hot gas is slower (~700 kms⁻¹) than in our fiducial models and has inner temperatures of about 10⁷ K, also in better agreement with the X-ray data. The resulting cloud kinematics is still in a relatively good agreement with the observed one, even though it is not as accurate as for our fiducial models (see Appendix B for more details). We stress that, to produce a similar wind, we need a large mass loading factor of about 1, which is higher than what is typically predicted by high-resolution simulations (e.g., Kim et al. 2020b; Schneider et al. 2020; Rathjen et al. 2021). η_M enters the model through the mass injection rate, defined as M_inj = n_MSFR, where we fix SFR = 0.1 M_⊙ yr⁻¹, consistent with observational estimates. However, the exact value of the SFR is still debated (see Section 5.2). If a higher SFR were assumed, a lower value of η_M would be required to reproduce hot wind properties consistent with X-ray observations. Overall, we conclude that our framework allows, taking into account all the observational and theoretical uncertainties, for a denser and cooler wind that would potentially be more in agreement with the X-ray constraints.

At Galactic latitudes higher than 1 degree, X-ray emission has recently been detected by eRosita (Predehl et al. 2021). The analysis of the eRosita spectra has shown that the X-ray emission is likely due to the combination of multiple components, some related to Galactic scales, like the eRosita bubbles and the MW hot circumgalactic medium (CGM; see Ponti et al. 2023; Locatelli et al. 2024) and others related to more local sources like the local hot bubble (Yeung et al. 2024). We note that to qualitatively better reproduce the emission measures of the Galactic X-ray components detected by eRosita, models with large mass loading factors (like those described above) are preferred to our fiducial best-fit models, which instead predict a less dense hot wind. However, given the uncertainties in both the observations and our model assumptions (see Sections 5.2 and 5.3), we leave a more thorough comparison with the eRosita observations for future work.

Fig. 9 Properties of the hot wind predicted by our best-fit models (in red; see main text for more details), compared to recent X-ray observational constraints from XMM-Newton (Ponti et al. 2019). From top to bottom, hot gas velocity, density, and temperature. The dashed curves show the predictions for a model where the hot mass loading factor ηM is forced to values close to 1 (see Section 5.1.1).

5.1.2 UV data

As already mentioned in Section 1, the cooler phases of the MW outflow have been detected not only through HI emission, but also extensively through UV absorption lines in the spectra of background quasars (QSOs) and halo stars (e.g., Fox et al. 2015; Bordoloi et al. 2017; Cashman et al. 2023). While the HI data are confined to latitudes lower than 10 degrees (but see Bordoloi et al. 2025), the sightlines where UV absorption was detected are all located at 10° < |b| < 50°, meaning that the two surveys trace two different areas of the MW wind, making it more challenging to assess whether the cold and warm phases of the outflow share the same physical origin.

Bordoloi et al. (2017) find that the velocities of the warm (T ~ 10^{4 5} K) clouds traced by the UV absorption go from about 250 kms⁻¹ at about 2.5 kpc from the disk down to 90 kms⁻¹ at a distance of 6.5 kpc. They model such kinematic signature with models of a decelerating outflow that is launched with initial intrinsic velocities of about 1000 kms⁻¹. Our findings, with cold clouds that are accelerated by a hot fast wind and reach radial intrinsic velocities of about 300 kms⁻¹ at 1-2 kpc (see Figure 8)¹⁰, therefore predict a distinct scenario compared to the outflow necessary to reproduce the warm clouds. This could mean that the two phases (warm and cold gas) might have a different kinematics, with that of the warm medium being potentially more similar to the kinematics of the hot wind. On the other hand, assuming that the clouds predicted by our models would survive up to heights larger than 2 kpc (but see Section 4.2) and would experience a deceleration when interacting with the preexisting material (see Section 5.2), our physical picture could potentially be able to reproduce at the same time both UV and 21-cm data. We plan to investigate this in a future study.

5.2 Main assumptions of the wind model

A fundamental assumption of this analysis is that the MW wind is powered by SN feedback. As already mentioned in Section 1, it is still debated whether the processes that originate the Fermi Bubbles and potentially also the wind entraining the cold HI clouds, are related to star formation or to the activity of Sgr A_* (e.g., Yang et al. 2022). The model of FB22, while developed for SN-driven winds, can in principle also be applied to an AGN-driven wind, by accordingly changing the injections of mass and energy in order to reproduce the expectations from the central SMBH. it is important to note that the model of FB22, originally based on the model of Chevalier & Clegg (1985), describes a steady-state wind, while AGN winds are blast wave-like (Sedov 1946; Taylor 1950) and are usually described by models that follow the time-dependent evolution of a shell of swept-up gas (see, e.g., King 2003; Faucher-Giguère & Quataert 2012; King & Pounds 2015; Richings & Faucher-Giguère 2018; Marconcini et al. 2025). While we acknowledge that this could be another channel to reproduce the observations of Section 2, developing such models for the MW wind is outside the scope of this work.

Throughout this paper, we have assumed that the wind is powered by the SN explosions within a spherical region at the center of the MW. Through Equations (13)-(16), we relate the mass and energy loading factors of the hot and the cold phase to the star formation rate. In our models, we assume SFR = 0.1 M_⊙/yr: this is consistent with observational estimates of the current star formation rate in the CMZ from both direct counting of young stellar objects and integrated light measurements (e.g., Barnes et al. 2017) and is expected to be relatively constant for at least the last 5 Myr and likely for a longer timescale (see Henshaw et al. 2023). While our assumption is therefore justified, variations of the adopted SFR by a factor of a few and a possible evolution throughout the 10 Myr analyzed in our work cannot be excluded (e.g., Nogueras-Lara et al. 2020, 2022). Such differences would most likely lead to slightly different values for the mass and energy loading factors that we inferred in Section 4. In addition to the value of the star formation rate, our model has the simplifying assumption of a spherical injection of mass and energy into the hot phase of the wind, while this region would be more realistically represented by a ring, resembling the structure of the CMZ (e.g., Armillotta et al. 2019; Sormani et al. 2020; Tress et al. 2020). A more sophisticated axisymmetric model for the hot gas phase (e.g., Sofue 2019; Nguyen & Thompson 2021, 2022) would naturally predict biconical instead of spherical outflows and different velocity, density and pressure profiles with respect to those found here and could be investigated in future studies. Our assumption is that the hot wind subtends a solid angle of 4π, while assuming a smaller angle would increase the wind densities, potentially leading to values closer to the X-ray observations (see Section 5.1.1). Previous toy models (Di Teodoro et al. 2018; Lockman et al. 2020) have shown that the kinematics of our population of HI high-velocity clouds can be described by a biconical wind with an opening angle of about 140°. This implies that the hot wind of our physical models needs to subtend at least such an angle in order to be responsible for the cloud entrainment. Reducing the opening angle to this value leads to an increase in normalization of the density (and pressure) profile of a factor 1.5, which would not change any of the findings of this paper.

Finally, we recall that in our model the hot wind is expanding in vacuum, while in reality the outflow is most likely interacting with the pre-existing hot gas components (e.g., Miller & Bregman 2016; Locatelli et al. 2024). This impact would likely modify the wind dynamics and properties. Sarkar (2024) roughly estimates the shock radius due to a galactic outflow expanding for a few megayears in the CGM, using the same expressions that are classically used to calculate the shock radius of a stellar wind on the interstellar medium (e.g., Weaver et al. 1977), to be about 2-3 kpc. Within this radius, one can expect the solution of the wind to be consistent with the vacuum-expanding one, which justifies our approach, given that all our findings are focused within roughly 2 kpc from the Galactic disk. Future works, trying to model the multiphase MW wind out to larger distances, should more thoroughly treat the impact between the outflow and the pre-existing CGM.

5.3 Uncertainties in the mixing between gas phases

One of the fundamental features of our model is the mixing of mass and momentum between the hot wind and the cold clouds due to turbulent mixing between them. An important approximation that we have adopted throughout this work is that the hot wind, even though it is continuously exchanging mass and momentum with the cold phase, maintains its properties fixed during the passage of the clouds and it is practically not affected by them. Here, we discuss how this simplification might influence our final findings.

We consider a radial multiphase model built entirely according to the full framework of FB22 (see the equations of their Figure 3) and assuming initial conditions that are consistent with our best-fit models of Section 4. Note that in here, the gravitational force influencing the cloud motion is only given by a logarithmic potential with circular velocity equal to 150 km s⁻¹, which represents the main difference with our models, where we instead consider a full axisymmetric representation of the MW potential (see Section 3.2). For this model, we used, for the hot phase, mass and energy loading factors equal to 0.15 and 0.5, respectively, and an injection radius equal to 150 pc; for the cold phase, we placed just outside the injection radius clouds with initial masses of 100 M_⊙ and initial kick velocities of 60 kms⁻¹, with a total cold gas mass loading factor equal to 0.2, and we assumed f_mix = 1/3. All the values above are in agreement with our fiducial best-fit models (see Figures 4 and 7). The hot wind properties without any exchange of mass and momentum from the clouds (as in the rest of this work) are shown by the gray curves in Figure 10. The solid purple curves show instead the properties of the hot gas once the exchange of mass and momentum is taken into account: one can see that the velocities decrease by less than a factor of two and they are still at about 1000 kms⁻¹, and that the profiles of density and pressure are flatter with respect to the previous case, in agreement also with the results of hydrodynamical simulations (e.g., Schneider et al. 2020). In the top panel, the dashed curve shows the predicted cloud radial velocity, showing that also in this case the clouds are significantly accelerated by the wind to velocities of 300 kms⁻¹ at a radius of 2 kpc. This is qualitatively in agreement (also considering the different gravitational potential of the two approaches) with the cloud velocities predicted by our fiducial models, implying that, despite the changes in the hot gas properties, the main findings of this work would not be altered.

As we have shown above and in Section 4, both the pressure exerted by the hot phase through the drag force and the accretion of high-velocity material are crucial to accelerating the HI clouds and reproducing the observed kinematics. In particular, the dimensionless parameter ξ determines whether the clouds experience an overall growth due to the gas mixing or whether they instead lose mass to the hot wind. We find that in our bestfit models ξ < 1 everywhere, which is consistent with the fact that our clouds are evaporating and becoming less massive with time (see Section 4). Another important parameter is f _mix, which determines the normalization of the mass loss/gain rate (Equation (11)) and that we left as a free parameter in our modeling. As we have seen in Section 3.2, this factor determines in practice the efficiency of the mixing between the cold and the hot phases of the outflow. Previous cloud-crushing simulations (e.g., Scannapieco & Brüggen 2015) have shown that clouds tend to be disrupted on timescales of N_cc cloud crushing times, where ${\dot{m}}_{\text{cl,loss}}\sim 1/N_{\mathrm{c}\mathrm{c}}{m}_{\mathrm{c}\mathrm{l}}/t_{\mathrm{c}\mathrm{c}}$. Comparing this with our Equation (10) implies that N_cc = 10/(3f_mix), so that for our best-fit models N_cc ~ 11. Interestingly, this value is roughly consistent with the results of previous works (e.g., Scannapieco & Brüggen 2015; Schneider & Robertson 2017), where N_cc ≳ 10 in the supersonic region (which is our region of interest, given that the cold clouds are injected outside of the sonic radius) and is instead lower in the subsonic region.

We stress that the above parametrization of N_cc, as well as the typical value of f_mix, is estimated from the results of previous hydrodynamical simulations. However, these simulations do not include, among others, physical effects such as cosmic rays (e.g., Ruszkowski & Pfrommer 2023; Armillotta et al. 2024) and thermal conduction (see, e.g., Brüggen & Scannapieco 2016; Armillotta et al. 2016; Afruni et al. 2023) and its potential interplay with magnetic fields (e.g., Kooij et al. 2021) nor variation in the hot background gas (e.g., Dutta et al. 2025): all these effects could potentially influence the amount of mixing between the different gas phases and therefore the cloud entrainment and survival. Our method of directly comparing analytical models with data can be used as a complementary way to infer the parameters (such as f_mix) of subgrid models (e.g., FB22; Huang et al. 2020) that can later be applied to cosmological simulations (e.g., Huang et al. 2022; Weinberger & Hernquist 2023; Smith et al. 2024) and can therefore be crucial for the future treatment of gas mixing in galaxy formation.

Fig. 10 Velocity (top), density (center), and pressure (bottom) profiles of a hot-wind model (solid curves) consistent with the best-fit models found in Section 4. In gray, we report a model where the exchange of mass and momentum with the clouds has no effect on the hot phase, while in purple we show the results from the same model, but taking into account the full framework of FB22. The dashed purple curve in the top panel shows the corresponding radial velocity of the cold clouds for such a model.

6 Summary and conclusions

In this work, we have created semi-analytical models of a multiphase Galactic wind, powered by SN feedback, with the aim of interpreting the observational data of McClure-Griffiths et al. (2013) and Di Teodoro et al. (2018), who through 21-cm emission found a population of HI high-velocity clouds close to the Galactic center. In our parametric models, which we built using GALPY (Bovy 2015) and the framework from Fielding & Bryan (2022), the hot gas phase tends to entrain and accelerate the cold clouds against the MW gravitational potential, through the effects of the drag force and of the exchange of mass (and therefore momentum) between the two phases. We infer the values of the model free parameters by performing a Bayesian analysis between the model outputs and the data.

Our main findings are the following:

1. We find a well constrained class of multiphase winds, driven by the star formation in the CMZ of the MW, that can successfully reproduce the locations, velocities and masses of the observed HI cloud population;

2. The hot gas phase of our best-fit models is well justified from a theoretical point of view, given that we infer mass and energy loading factors in very good agreement with expectations from recent high-resolution hydrodynamical simulations. Our predicted hot winds are fainter and hotter compared to current X-ray constraints: this could either mean that we are modeling a different component to that detected in the observations, or that larger mass injection rates are needed;

3. The cold clouds have masses and (roughly) radii in agreement with the observational estimates and are accelerated by the interactions with the hot wind to velocities of 300400 kms⁻¹. We also find that, in order to reproduce the observations, the clouds need to overall lose mass to the hot component, with the least massive clouds having lost more than 70% of their initial mass after 10 Myr.

This work therefore provides a physically motivated and self-consistent model that is able to reproduce the kinematics of this population of HI high-velocity clouds. With this analysis, we have found that SN feedback (even though we cannot exclude the influence of the central AGN) could be a natural driver of the multiphase outflow at the center of our Galaxy, and it is therefore a crucial process to take into account for the creation of the Fermi Bubbles.

Acknowledgements

We thank Kartic Sarkar, Gabriele Ponti, Rongmon Bordoloi, Drummond Fielding, Alankar Dutta, Gabriele Pezzulli and Filippo Fraternali for helpful discussions that have contributed in developing some of the ideas presented in this work. AA and EDT were supported by the European Research Council (ERC) under grant agreement no. 101040751. LA was supported by the INAF Astrophysical fellowship initiative.

References

Appendix A Ballistic model

In this appendix, we investigate whether a ballistic model would be able to reproduce the kinematic data described in Section 2 and shown in Figure 1, to assess whether or not a hot entraining wind is needed to explain the observed kinematic pattern. For this purpose, we modeled the cloud orbits using GALPY as in Section 3.2, but including only the MW gravitational potential and no other additional force due to the influence of the hot wind. We assumed that all the cold clouds start within the injection radius (r_inj, which in this case obviously refers only to the cold gas), with initial (R_in, z_in) for each orbit chosen exactly as in Section 3.2. As for the initial cloud velocities, we assumed an initial kick velocity, υ_kick, as in our fiducial models. Differently from Section 3.2, we did not assume that υ_kick is entirely vertical (that was previously justified by the fact that the influence of the hot wind would dominate the cloud dynamics), but we assumed an angle, θ_max, that determines the extent of the biconi-cal outflow aperture. For each orbit, we indeed selected a random θ uniformly sampled between 0 and θ_max. We additionally chose for each orbit a random angle ϕ between 0 and 2π, such that (see Afruni et al. 2021): $$\begin{array}{r}{v}_{\mathrm{k}\mathrm{i}\mathrm{c}\mathrm{k},\mathrm{R}}=v_{\mathrm{k}\mathrm{i}\mathrm{c}\mathrm{k}}\sin \theta \cos \varphi ,\\ v_{\mathrm{k}\mathrm{i}\mathrm{c}\mathrm{k},\mathrm{T}}=v_{\mathrm{c}\mathrm{i}\mathrm{r}\mathrm{c}}+v_{\mathrm{k}\mathrm{i}\mathrm{c}\mathrm{k}}\sin \theta \sin \varphi ,\\ v_{\mathrm{k}\mathrm{i}\mathrm{c}\mathrm{k},\mathrm{z}}=v_{\mathrm{k}\mathrm{i}\mathrm{c}\mathrm{k}}\cos \theta .\end{array}$$(A.1)

In this framework, the free parameters of our models are therefore r_inj, θ_max and υ_kick. For each choice of these free parameters we modeled 20 different orbits for 10 Myr, with initial conditions as explained above, and we compared the model clouds populating these orbits (created exactly in the same way as in Section 3.2.2) with the observed HI components, using the likelihood defined by equation (12). In particular, we ran a Bayesian analysis with DYNESTY, with the following priors: a flat prior on the injection radius 0.05 < r_inj/kpc < 0.3 (in order to be consistent with the CMZ extent, see also Section 3.3); a flat prior in θ_max between 0 (vertical outflow) and π/2 (isotropic outflow); a flat prior in the kick velocity 0.5 < υ_kick,₁₀₀ < 5.

We show the results of this analysis in Figure A.1. In the top panel, we report the corner plot of the posterior distributions of the three free parameters: one can see that in order to reproduce the data the models need large injection radii (r_inj ≳ 250 pc), relatively large outflow apertures (θ_max ~ π/3) and large initial kick velocities (υ_kick ~ 330 kms⁻¹). In particular, both the outflow velocity and aperture are very well constrained and exhibit narrow posterior distributions. In the bottom panel of Figure A.1 we show instead how these best-fit models perform in reproducing the data. Here, we report the cloud velocity as a function of the Galactic latitude, with the data in orange and the bestfit models in teal (similar to Figure 6). One can immediately see that, contrary to our fiducial multiphase models, a ballistic model is not able to reproduce the observed kinematic pattern of acceleration. The models show indeed, as expected, the opposite behavior of deceleration, given that the gravitational force tends to slow down the clouds that are initially ejected with very high velocities. This comparison therefore demonstrates how the motion of the HI clouds detected in the vicinity of the Galactic center cannot be explained only by gravity.

Note that taking into account the impact of the cold outflow with a static or slowly rotating hot CGM (e.g., Miller & Bregman 2015; Hodges-Kluck et al. 2016; Locatelli et al. 2024, see also Section 5.2), similar to the models implemented in Afruni et al. (2021, 2022), would make the model outputs even more inconsistent with the data, as it would increase the amount of deceleration of the clouds. We therefore conclude that a fast-moving hot wind is needed to reproduce the observed kinematic pattern of the HI high-velocity clouds.

Fig. A.1 Main results of the ballistic models presented in Appendix A. Top: corner plot showing the posterior distributions of the three free parameters; the vertical lines mark the positions of the median values of the posteriors and of the 2 - σ uncertainties. Bottom: comparison between the outputs of the best-fit models (teal) and the observations (orange).

Appendix B Massive hot wind model

In Section 5.1 we have discussed a model where we imposed the hot wind mass loading factor to be equal to 1, by adopting in our Bayesian analysis a narrow Gaussian prior centered in log η_M = 0. This leads to a slower, denser, and cooler wind, in better agreement with the observational constraints from the X-ray chimneys data (see Figure 9). Even though our preferred models have lower mass loading factors, more consistent with theoretical predictions (see discussion in Section 5.1.1), we show here that if we considered a more massive wind model the observed HI cloud kinematics would still be reproduced by our framework.

Fig. B.1 Corner plot showing the posterior distributions of the five free parameters of the model presented in Section 5.1.1 and Appendix B, where a narrow Gaussian prior centered in log ηM = 0 is adopted for the hot wind mass loading factor. As for our fiducial model (see Section 3.3), we also adopted a Gaussian prior for the energy loading factor, centered in log ηE = −0.4. The vertical lines mark the positions of the 2.5, 50 (solid line) and 97.5 percentiles of the 1-dimensional posterior distributions.

In the two panels of Figure B.2, we show the outputs of a model where the values of the 5 parameters are fixed to the median values of the posterior distributions obtained with the Bayesian analysis (Figure B.1). In the top panel, we show the evolution of the total radial velocities (obtained by summing in quadrature the three velocity components) of 20 typical clouds in this model. Note how, as for our fiducial models, the clouds are accelerated by the interactions with the hot wind, with the lower mass clouds reaching higher velocities compared to the more massive ones. Compared to the scenario presented in Section 4, where the clouds keep accelerating with increasing distance from the center, here they appear to reach a maximum velocity of about 300 km s⁻¹ (for the lower mass clouds). This effect is likely due to the lower velocities of the hot wind in this model, as can be seen in Figure 9. As a result (see bottom panel of Figure B.2), this model does not reproduce the few clouds at high latitudes with velocities approaching 400 km s⁻¹, which were instead recovered by our fiducial model. However, except for this small discrepancy, the cloud kinematic distribution is still consistent with the observed one. This leads us to conclude that a denser wind would still be able to overall explain the kinematics of this population of HI high-velocity clouds close to the Galactic center.

Fig. B.2 Outputs of the best-fit model assuming a massive hot wind (see Figure B.1 and Appendix B). Top: Radial velocities of 20 typical clouds as a function of the distance from the center, color-coded by the cloud initial mass. Bottom: Comparison between the outputs of the bestfit models (teal) and the observations (orange).

All Tables

All Figures

Fig. 1 Population of HI high-velocity clouds from the surveys of McClure-Griffiths et al. (2013) and Di Teodoro et al. (2018), taken with the ATCA and the GBT. Left: galactic latitude-longitude map. The data-points show the positions of the clouds and are color-coded by their vLSR. The white mask shows the region surveyed by the ATCA and the GBT. Right: HI cloud latitudes as a function of their velocities. The kinematic pattern shows a signature of acceleration (see main text and Lockman et al. 2020). The points are color-coded by the cloud mass and the full 1D mass distribution is shown in the inset panel.

In the text

Fig. 2 Behavior of the multiphase wind for three different choices of models (see main text): model A (top row), model B (central row), and model C (bottom row). On the left, hot gas velocities (solid curves) and densities (dashed curves) are shown as a function of the distance, r. The right panels show instead the cold cloud orbits (solid curves), indicating in particular the effects of drag (dashed curves) and mixing (dash-dotted curves) on the original ballistic orbit (dotted curves). The symbols mark the positions of the modeled cloud at the beginning (circles), after 5 Myr (squares) and after 10 Myr (rhomboids).

In the text

Fig. 3 Cloud mass evolution as a function of the distance from the center, for models A, B, and C (see Figure 2). In all cases, the cloud starts with a mass of 5 × 102 M⊙, but in model A it continuously loses its mass and slowly evaporates into the background medium due to hydrodynamical interactions with the hot wind, while in the other two models it initially increases its mass due to the cooling of the mixing layers.

In the text

	Fig. 4 Corner plot showing the posterior distributions of the five free parameters of our models. The vertical lines mark the positions of the 2.5, 50 (solid line) and 97.5 percentiles of the one-dimensional posterior distributions.
In the text

	Fig. 5 Comparison of the outputs of the best-fit models with the observational data. The cloud kinematics is shown as a function of Galactic longitude and latitude, with models on the top and data on the bottom. The one-dimensional distributions of the observed (orange, solid line) and the model (teal, dashed line) longitudes are reported on top.
In the text

	Fig. 6 Comparison of the outputs of the best-fit models with the observational data. The model (teal) and data (orange) line-of-sight velocities are shown as a function of the Galactic latitude. The marginal onedimensional distributions are plotted along the two axes for both the data and the model.
In the text

	Fig. 7 Mass (top) and energy (bottom) loading factors as a function of the galactocentric radius, r, for a class of best-fit models (see Section 4.2 for more details), for the hot gas (red) and the cold clouds (blue). The vertical dashed line marks the location of the median injection radius.
In the text

	Fig. 8 Intrinsic properties of the HI clouds predicted by our best-fit models. Top, from left to right: cloud orbits, cloud masses as a function of r, and cloud radii as a function of r. Bottom: three components of the cloud velocities. The different curves are colored based on the initial cloud mass, from lower masses in cyan to higher masses in magenta.
In the text

Fig. 9 Properties of the hot wind predicted by our best-fit models (in red; see main text for more details), compared to recent X-ray observational constraints from XMM-Newton (Ponti et al. 2019). From top to bottom, hot gas velocity, density, and temperature. The dashed curves show the predictions for a model where the hot mass loading factor ηM is forced to values close to 1 (see Section 5.1.1).

In the text

Fig. 10 Velocity (top), density (center), and pressure (bottom) profiles of a hot-wind model (solid curves) consistent with the best-fit models found in Section 4. In gray, we report a model where the exchange of mass and momentum with the clouds has no effect on the hot phase, while in purple we show the results from the same model, but taking into account the full framework of FB22. The dashed purple curve in the top panel shows the corresponding radial velocity of the cold clouds for such a model.

In the text

	Fig. A.1 Main results of the ballistic models presented in Appendix A. Top: corner plot showing the posterior distributions of the three free parameters; the vertical lines mark the positions of the median values of the posteriors and of the 2 - σ uncertainties. Bottom: comparison between the outputs of the best-fit models (teal) and the observations (orange).
In the text

Fig. B.1 Corner plot showing the posterior distributions of the five free parameters of the model presented in Section 5.1.1 and Appendix B, where a narrow Gaussian prior centered in log ηM = 0 is adopted for the hot wind mass loading factor. As for our fiducial model (see Section 3.3), we also adopted a Gaussian prior for the energy loading factor, centered in log ηE = −0.4. The vertical lines mark the positions of the 2.5, 50 (solid line) and 97.5 percentiles of the 1-dimensional posterior distributions.

In the text

	Fig. B.2 Outputs of the best-fit model assuming a massive hot wind (see Figure B.1 and Appendix B). Top: Radial velocities of 20 typical clouds as a function of the distance from the center, color-coded by the cloud initial mass. Bottom: Comparison between the outputs of the bestfit models (teal) and the observations (orange).
In the text