Resolving the Unresolved Galactic Winds in Multiphase Models. I. Methodology and Application

Our sample includes 22 local star-forming galaxies jointly compiled from the Cosmic Origins Spectrograph (COS) Legacy Archive Spectroscopy Survey (CLASSY) atlas (D. A. Berg et al. 2022b; B. L. James et al. 2022) and the sample presented by T. M. Heckman et al. (2015). These galaxies were observed with the Hubble Space Telescope (HST) COS, forming the largest and highest-quality UV spectroscopic dataset of local star-forming galaxies. These objects span a wide range of important galaxy properties, including stellar mass (log(M_⋆) ∼ 6–11 M_⊙), star-formation rate (SFR; ∼0.01–100 M_⊙ yr⁻¹), metallicity (12 + log(O/H) ∼ 7–9), and [S ii]-based electron density (n_e ∼ 10¹–10³ cm⁻³). The COS observations, given in an aperture¹⁷ of 25, result in spatially integrated spectra covering the far-UV wavelength range from ∼1200 to 2000 Å. To improve the signal-to-noise ratio, the spectra are resampled into bins of 0.18 Å, corresponding to a resolving power of R ∼ 6000–10,000, or velocity resolution of 50–30 km s⁻¹ from the blue to red ends. Further details on the program design and data reduction can be found in previous publications (D. A. Berg et al. 2022b; B. L. James et al. 2022; T. M. Heckman et al. 2015). The CLASSY data were obtained from the Mikulski Archive for Space Telescopes (D. Berg et al. 2022a).

Given the dedicated UV dataset, cool outflow properties in these galaxies have been studied in various publications (e.g., X. Xu et al. 2022, 2023b; W. Hu et al. 2023; M. Huberty et al. 2024; K. S. Parker et al. 2024). We take the one that focus on the rest-frame UV absorption lines (X. Xu et al. 2022). In brief, outflow signatures are detected in lines including O i λ1302, C ii λ1334, and the Si ii multiplet (λλλλλ1190, 1193, 1260, 1304, 1526), and from higher-ionization transitions, e.g., Si iii λ1206, and Si iv λλ1393, 1402 (see examples in the left panel of Figure 1). Then the velocity-dependent ionic column density (dN_ion/dv) is measured using the partial coverage method (PCM, e.g., D. S. Rupke et al. 2005; C. L. Martin & N. Bouché 2009; J. Chisholm et al. 2016; T. E. Rivera-Thorsen et al. 2017). These column densities are then fed into the Cloudy photoionization models (v.17.01; G. J. Ferland et al. 2017) to constrain the total hydrogen column density per velocity (dN_H/dv). An example of the dN_H/dv profile is shown as the blue lines in the right panel of Figure 1. After that, X. Xu et al. (2022) adopt spherical outflow models to estimate the average mass, momentum, and energy outflow rates (, , and , respectively), which represent the spatially integrated values within the COS aperture. Finally, the cool outflows’ mass-loading factors (η_M,cool) are calculated as /SFR, where SFR is derived from spectral energy distribution (SED) fittings described in D. A. Berg et al. (2022b).

Zoom In Zoom Out Reset image size

To model galactic winds, one must capture both their multiphase and multiscale character—that is, how distinct temperature phases interact as the flows propagate from the galaxy into the CGM. In addition, to compare with observational results (Section 2.1), the model must remain flexible enough to span diverse galaxy properties, including SFR, M_⋆, and the outflow mass-loading factors in different phases (η).

With these requirements in mind, we adopt the state-of-art analytic wind model developed by D. B. Fielding & G. L. Bryan (2022, hereafter FB22). In this framework, the galactic wind consists of a hot, volume-filling component (the hot wind) and a cool component composed of embedded outflowing clouds (the cool outflows). As these two phases coevolve, they exchange mass, momentum, and energy through turbulent radiative mixing layers (TRMLs), a process motivated by results from high-resolution simulations (e.g., S. Ji et al. 2019; M. Gronke & S. P. Oh 2020; N. Mandelker et al. 2020; M. W. Abruzzo et al. 2024; T. Marin-Gilabert et al. 2025). The fate of the cool clouds is governed by the relative timescales of the radiative cooling (t_cool) in TRMLs surrounding each cloud and their mixing with the hot wind (t_mix = R_cl/V_turb), where V_turb is the turbulent velocity. Large clouds with t_cool < t_mix are able to condense additional mass from the hot wind and grow. In contrast, smaller clouds with t_cool > t_mix are shredded by turbulent mixing with the hot wind. This ongoing mass and energy exchange not only determines the survival of the cool clouds, but also modifies the global structure of the hot wind relative to classical single-phase models.

The model considers that mass and energy are injected uniformly inside the star-forming region radius (r < r_*), and cooling and gravity are not included within this region. Beyond this radius, radiative cooling and gravity are included, which shape the coupled evolution of the wind and clouds as they traverse the surrounding medium. Cooling curves from R. P. C. Wiersma et al. (2009) are adopted, including heating from the z = 0 UV background. The flow is treated as one dimensional, with forces confined to the radial direction. In addition, we assume a spherically symmetric galactic wind in the model, which is compatible with the large outflow solid angle (Ω) inferred from our observations (X. Xu et al. 2022). The temperature of the cool outflow clouds is set to be T = 10⁴ K, which also matches the one taken for the observed UV absorption lines. We assume the hot wind has a metallicity of 2 Z_⊙ (e.g., L. A. Lopez et al. 2020), while we adopt the observed gas-phase metallicity for the cool outflows (D. A. Berg et al. 2022b).

In what follows, we show how the multiphase, multiscale FB22 framework is applied to fit single-phase, spatially integrated outflow measurements. The first task is to select observables that capture the properties of the cool outflow clouds and, via these, indirectly constrain the hot wind and the outflow’s radial structure.

To characterize the cool outflow clouds, the key observables are their mass, momentum, and energy outflow rates, which track the strength and impact of feedback. We start with the observational analysis of outflows. Under a spherical thin-shell outflow model, mean rates follow directly from the velocity-resolved column density profile (e.g., X. Xu et al. 2022):

where μ is the mean atomic mass per proton (∼1.4), m_p is the proton mass, R₀ is the mean outflow radius from which the outflow begins, v is the mean outflow velocity, and Ω is the global covering factor, which approaches 4π as described in X. Xu et al. (2022). We adopt R₀ = 2 × r₅₀, where r₅₀ is the UV half-light radius of the star-forming region. This choice is consistent with previous studies (e.g., P. L. Shopbell & J. Bland-Hawthorn 1998; T. M. Heckman et al. 2015; X. Xu et al. 2023b) and corresponds to an outflow that is launched at a radius enclosing ∼90% of the starburst. The remaining free parameter is the velocity-dependent total hydrogen column density of the clouds (N_H(v)).

Given the moderate spectral resolution of HST/COS (30–50 km s⁻¹ in our case), our galaxies enable velocity-resolved measurements of N_H(v). An example of it appears as the blue line in the right panel of Figure 1. We fit three summary quantities of N_H(v) that capture its normalization, characteristic velocity, and width of the profile. This choice balances observability and fidelity: each quantity is straightforward to measure and together they reproduce the profile’s overall shape. Specifically, we use (1) the integrated hydrogen column of the outflowing clouds (N_H,cl), obtained by integrating over v < 0 km s⁻¹; (2) the bulk velocity, v_cl, defined as the velocity at which dN_H/dv reaches its maximum; and (3) the half-width at half maximum on the blueshifted side (HWHM_cl), i.e., the width of the high-speed wing with v < v_cl. We do not use the full width because the FB22 model produces only outflows, so the modeled profiles occupy only the blueshifted side (v < 0).

The FB22 model directly yields the velocity and number density of outflow clouds at each radius (v_cl(r) and n_H,cl(r)). From that, we have

We then measure N_H,cl, v_cl, and HWHM_cl from this modeled dN_H/dv profile and fit them to the ones measured from observations.

Next, we identify the key free parameters varied in the model. FB22 have identified the primary ones that govern the outputs. On the one hand, these include galaxy-scale properties, e.g., SFR, the depth of the gravitational potential (quantified by M_⋆ or v_cir), and the size of the star-forming region (r_* ∼ r₅₀ measured in the UV). For our sample, these quantitatives are measured directly from observations (D. A. Berg et al. 2022b; X. Xu et al. 2022); therefore we fix them to the observed values (Table 1). On the other hand, there are parameters that mark the initial status of the winds and outflows (i.e., their values at r_*). These include the initial hot-phase mass-loading factor (η_M,hot,0 ≔/SFR), the initial cool-phase mass-loading factor (η_M,cool,0 ≔ /SFR), and the initial outflow cloud mass (M_cl,0). We select them because previous analyses have shown that UV absorption lines place tight constraints on the characteristic η_M,cool and M_cl (X. Xu et al. 2022, 2023b). Additionally, although the hot phase is not directly detected, its interaction with the cool phase may allow us to set limits on η_M,hot,0.

Table 1. Model Parameters and Results for the Sample

Object	SFR	r50	log(M⋆)	vcl	HWHMcl	log(NH,cl)	log(ηM,hot,0)	log(ηM,cool,0)	log(Mcl,0)	χ2
	(M⊙ yr−1)	(kpc)	(M⊙)	(km s−1)	(km s−1)	(cm−2)			(M⊙)
(1)	(2)	(3)	(4)	(5)	(6)	(7)	(8)	(9)	(10)	(11)
J0021+0052	11.75	0.45	9.09	366	158	20.71	0.36	–0.29	6.27	3.38
J0036–3333	10.23	0.12	9.14	150	83	20.45	–0.30	–1.00	6.09	0.36
J0055–0021	42.66	0.57	9.56	466	219	20.92	0.18	–0.48	6.23	9.46
J0150+1308	31.62	0.63	9.80	175	125	20.63	–0.53	–0.50	5.99	0.80
J0823+2806	30.20	0.25	9.38	199	164	20.44	–0.20	–0.89	6.02	3.36
J0938+5428	11.22	0.51	9.15	153	36	20.78	–0.08	–0.01	6.18	0.27
J0942+3547	0.17	0.10	7.56	100	58	20.03	0.44	0.12	4.62	0.61
J1024+0524	1.62	0.27	7.88	103	56	20.91	0.15	0.42	6.30	0.27
J1025+3622	10.96	0.79	8.87	149	101	20.79	–0.20	0.11	6.29	0.09
J1112+5503	39.81	0.47	9.59	389	126	20.46	–0.35	–0.59	3.75	0.69
J1144+4012	32.36	0.91	9.89	177	118	20.64	–0.40	–0.32	6.13	0.04
J1150+1501	0.05	0.07	6.84	71	32	20.54	0.68	0.82	5.17	0.75
J1157+3220	9.33	0.65	9.04	161	47	20.13	0.24	–0.60	6.04	0.04
J1200+1343	5.62	0.23	8.12	300	198	20.88	0.44	–0.13	6.19	2.16
J1253–0312	3.63	0.39	7.65	123	39	20.85	–0.12	0.24	6.12	1.19
J1414+0540	5.00	0.36	9.69	178	107	20.39	–0.56	–1.02	6.26	1.26
J1416+1223	37.15	0.27	9.59	382	111	20.07	–0.35	–0.87	3.52	0.02
J1428+1653	16.60	1.04	9.56	136	81	20.61	–0.59	–0.09	6.29	0.75
J1429+0643	26.30	0.44	8.79	238	151	20.74	–0.50	–0.33	5.04	0.78
J1448–0110	2.45	0.12	7.61	132	84	20.36	0.14	–0.50	5.66	1.47
J1521+0759	8.91	0.48	9.00	139	67	20.33	–0.36	–0.31	6.23	0.69
J2103–0728	19.50	0.24	10.28	200	72	20.08	–0.03	–1.05	5.43	0.04

Note. We show galaxy properties in columns (2)–(4) and the galactic-wind and outflow properties in columns (5)–(11). Among them, columns (2)–(7) are fixed and adopted to constrain the model parameters in columns (8)–(10). The galaxy properties are measured from HST/COS spectra and images as described in X. Xu et al. (2022). The initial parameters (i.e., at r = r₅₀) for the wind and outflow properties from the best-fit FB22 model (see Section 2.3). We vary three parameters (columns (8)–(10)) to fit the observed three moments of the dN_H/dv profiles (columns (5)–(7)). The 1σ error bars are shown with the values, which are derived from the Markov Chain Monte Carlo (MCMC) chain. The final column shows the reduced χ² values of the fit.

Download table as:  ASCIITypeset image

Overall, we fit the observed v_cl, HWHM_cl, and N_H,cl for each galaxy by varying three free parameters in the FB22 model, specifically with log(η_M,hot,0) ∈ (−1.5, 1.5), log(η_M,cool,0) ∈ (−1.5, 1.5), and log(M_cl,0/M_⊙) ∈ (1, 6.5). We also limit M_cl,0 < 0.1 M_⋆ for each galaxy. The observed galaxy and outflow properties are listed in the first seven columns of Table 1.

We perform the fits using a Bayesian framework implemented via the emcee MCMC sampler (D. Foreman-Mackey et al. 2013). We sample the posterior probability distribution using N = 72 walkers for a total chain length of , where is the largest integrated autocorrelation time across all parameters. This is to ensure adequate convergence and independent sampling. The resulting posterior distributions provide both median estimates and 16th–84th percentile credible intervals for all parameters. The fitting results are listed in the last four columns of Table 1 and discussed in Sections 3 and 4.

While the FB22 model represents a significant step forward in developing simple and physically motivated multiphase wind models, it necessarily assumes time independence and a specific form for the cloud–wind interaction, both of which are approximations to the full problem. In addition, several parameters in that model are held fixed at fiducial values because the present data do not have the leverage to constrain the full parameter space. These include the normalization of the turbulent-mixing-driven cloud-mass growth rate () and the parameter setting the strength of turbulence in the TRML relative to the cloud–wind relative velocity (f_turb = 0.1). Therefore, the results presented in the following sections should be interpreted as effective, model-conditional constraints. Future papers will relax some of these assumptions and explore a broader and more flexible model hierarchy to assess how robust these inferences are to changes in the underlying physics.

First, we assess how well the three free parameters (η_M,hot,0, η_M,cool,0, and M_cl,0) are constrained by the data. We explore it with a three-dimensional grid of models. Figure 2 visualizes the parameter space: each panel scans two parameters on the axes, while the third is fixed to the galaxy’s best-fit value. The first three rows show model predictions for v_cl, HWHM_cl, and N_H,cl as the color scales, indicating the constraining power of each observable. The bottom row shows the reduced χ² surfaces. Black stars mark the best-fit model location for J0021+0052. Blank regions denote parameter combinations that produce nonphysical solutions with the FB22 model and are excluded from further consideration.

Zoom In Zoom Out Reset image size

We note that η_M,cool,0 is tightly constrained, with uncertainties < 0.1 dex. This is most evident in panel (10) of Figure 2. Such tight constraints are expected because the fitted moments of dN_H/dv can be mapped directly to (Equation (1)) and, for a fixed SFR, to η_M,cool. By contrast, M_cl,0 is less tightly—but still moderately—constrained, with typical uncertainties of 0.2–0.3 dex (e.g., panels (3), (6), and (12)). This weaker constraint arises because M_cl scales with N_H but also depends on additional factors, including the outflow hydrogen number density (n_H) and and how outflow clouds shield each other in the line of sight (LOS; see Section 4.4 in X. Xu et al. 2023b).

η_M,hot,0 is moderately constrained with uncertainties of ∼0.3 dex (e.g., panels (1), (4), (10) of Figure 2). Among the three observables, v_cl carries more leverage when η_M,hot,0 > 1, whereas HWHM_cl is more informative when η_M,hot,0 < 1. These can be seen from the changes of color gradients in the horizontal directions in panels (1) and (4), respectively. Overall, this supports the expectation that, although the observations probe only the cool outflow clouds, the hot winds can be constrained indirectly because the two phases interact as they propagate out of the galaxy. We reach a similar conclusion for other galaxies in our sample, which we discuss more in Section 4.

We note that Ω = 4π is assumed in our models, which is a statistical conclusion from X. Xu et al. (2022) for local starburst galaxies. Since and the corresponding η_M are linearly proportional to Ω (see Equation (1)), our derived η_M,cool,0 and η_M,hot,0 should be viewed as angle-averaged quantities. For example, if a galaxy actually have Ω = 2π, then both η_M,cool,0 and η_M,hot,0 should drop by a factor of 2. We also find the model fit with a smaller Ω yield a larger . This is likely because reducing Ω decreases the total number of clouds intersecting the LOS, which can be compensated by larger cloud masses to reproduce the observed column density distribution.

Having established how the parameters are constrained, we now ask whether the spatially integrated spectra can recover spatial information about the outflows. This is plausible because dN_H/dv derived from absorption lines should already encode spatial information about the outflow structure (e.g., C. Carr et al. 2018; M. Huberty et al. 2024; Z. Peng et al. 2025; C. A. Carr et al. 2025; X. Xu et al. 2025). Specifically, clouds at smaller radii are generally denser and move more slowly, contributing absorption at lower velocities. Conversely, clouds at larger radii are more diffuse and travel faster, producing absorption at higher velocities. High spectral resolution in absorption-line profiles therefore helps constrain the radial distribution of the outflowing gas. We illustrate this mapping (from velocity to radius) by annotating the model-predicted radius along the top axis of the right panel of Figure 1.

We note that the velocity–radius mapping is strictly one to one only when the bulk velocity field is monotonic. In more general cases with nonmonotonic flows, the inversion from velocity to radius is no longer unique, as gas at multiple radii can contribute to the same observed velocity. Nevertheless, the absorption profile still encodes useful information about the global distribution of the outflowing gas (e.g., the total column density and characteristic velocity scales), although its association with specific radii becomes model dependent. For the galaxies modeled in this paper, all systems are well fit with monotonic velocity distributions (Section 4.2), so the mapping remains well defined. We also emphasize that the FB22 framework does not require monotonic velocity profiles; monotonicity simply emerges from the best-fit solutions for this sample.

To clarify the interaction between the hot and cool phases and their radial structure, we present radial profiles of key quantities in Figure 3. The blue curves come directly from the FB22 model using the best-fit solution for J0021+0052, where solid and dashed lines denote the hot wind and the cool outflow clouds, respectively. Panel (1) shows a rapid acceleration of the cool outflow from r₅₀ (0.45 kpc for J0021+0052) to ∼1 kpc, followed by a plateau at larger radii. The hot-wind velocity, v_wind, also exhibits a slight deceleration beyond 1 kpc. These trends are consistent with continued mass entrainment from the hot phase by the cool clouds, leading to an increase in the cloud mass (panel (4)). Notably, the hot wind remains supersonic at all radii. This behavior matches the criterion that D. D. Nguyen et al. (2024) suggest for one-dimensional time-steady analytic wind solutions to reproduce three-dimensional time-dependent calculations. In red, we overlay the observed outflow velocity, v_cl (see Section 2.3), measured from the spatially integrated HST/COS spectrum (X. Xu et al. 2022), with shaded bands indicating the 1σ uncertainties. We find it intersects the modeled v_cl(r) at ∼0.7–0.8 kpc, close to 2 r₅₀ for this galaxy.

Zoom In Zoom Out Reset image size

For the whole sample, we perform similar comparisons (detailed in Section 4) and find that the intersections consistently occur between 1 and 2 r₅₀. This suggests that the spatially integrated velocities derived from our observations correspond to the resolved velocities at a characteristic radius. This range is consistent with previous work that found outflow sizes extending to 1–2 r₅₀, including S. F. Newman et al. (2012), who directly measured outflow sizes from optical emission lines, and X. Xu et al. (2023b) who indirectly derived outflow sizes using UV absorption lines.

The bottom row shows the radial profiles of the mass, momentum, and energy outflow rates. A clear transfer of mass, momentum, and energy from the hot phase (solid blue) to the cool phase (dashed blue) can be seen. Taken together, interaction between the cool and hot phases significantly shapes the properties of both. The velocity–radius mapping also enables us to recover spatial information from integrated spectra. In the next section, we extend this analysis to the full sample.

For all 22 galaxies in our sample, their observed galaxy and outflow measurements are listed in Table 1, where the best-fit parameters (η_M,hot,0, η_M,cool,0, and M_cl,0) and reduced χ² values are summarized in the last four columns.

Overall, we find most galaxies fit well with the models (χ² < 1.5 for 18 out of 22 galaxies), while the remaining four have a relatively worse fit, which we discuss in Section 5.2. We show the distributions of the fitting parameters in Figure 4 as histograms. We also overlay posterior probability density functions (PDFs) from the MCMC to reflect parameter uncertainties. In the left panel, we find ∼60% of our galaxies prefer a hot-phase initial mass-loading factor (η_M,hot,0) < 1.0. This is consistent with previous theoretical studies of η_M,hot, e.g., in D. Zhang et al. (2014). They highlight that heavily mass-loaded winds cannot be described by the adiabatic R. A. Chevalier & A. W. Clegg (1985, hereafter CC85) model because they become strongly radiative. They suggest η_M,hot ≲ 1 when SFR ≳ 10 M_⊙ yr⁻¹, which is mostly true for our galaxies matching this criterion. In addition, there are three galaxies with a large η_M,hot,0 ≳ 3 (J0942+3547 and J1150+1501, and J1200+1343). This is likely due to their very low SFRs (<0.2 M_⊙ yr⁻¹) derived from SED fitting (D. A. Berg et al. 2022b), which can have larger relative uncertainties. Though this assumption was found to provide satisfying fits to the broadband photometry of our galaxies, more complex star formation histories (e.g., bursty SF) may explain why we observe outflow signatures in these galaxies with a very low SFR.

Zoom In Zoom Out Reset image size

For η_M,hot, V. Pandya et al. (2021) found the following relationship for galaxies from the FIRE-2 simulations based on spatially integrated measurements:

This suggests that η_M,hot ≲ 3 for the mass range of our galaxies (i.e., M_⋆ ≳ 10⁷ M_⊙), which is also consistent with what we find in Table 1.

In the middle panel of Figure 4, we show the cool-phase initial mass-loading factors (η_M,cool,0). We find most of our galaxies (∼80%) prefer η_M,cool,0 < 1, with a PDF that peaks at ∼0.1. This is expected since the cool outflow clouds initially are at a relatively low velocity and mass when they are not accelerated by the hot wind.

In the right panel, we show that the initial mass of the outflow clouds (M_cl,0) favors high values approaching 10^6.5 M_⊙. We caution that the current model assumes that all outflow clouds have the same mass, which is likely unrealistic. Therefore, we leave the interpretation of M_cl to the subsequent paper, when we will update the model with more physically motivated mass distributions of cool outflow clouds.

As in Section 3.2, we present our best-fit model-predicted radial distributions of the outflow properties for all galaxies in Figures 5, 6, and 7. Each curve corresponds to one galaxy and is color coded by the M_⋆ of that galaxy. The horizontal axis shows the radius normalized by its r₅₀.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

In Figure 5, one empirical finding is that the cool-phase outflow velocity, v_cl, increases monotonically with radius for all galaxies (solid lines). This behavior is consistent with the cool outflow’s acceleration driven by the hot wind through mainly momentum transfer. The nondetection of decelerating outflows is likely because such outflowing clouds lose mass and persist for shorter timescales, making them harder to be observed. There is a strong trend in which more massive galaxies show a higher v_cl, which is known through the outflow scaling relationships in previous studies (e.g., T. M. Heckman et al. 2015; X. Xu et al. 2022).

We also overlay the observed spatially resolved velocity profile from M82 as a black line with its 1σ error in gray shaded regions. This profile approximates a β law (X. Xu et al. 2023a):

where v_∞ is the terminal outflow velocity and β is the power-law index. We find that the M82 profile exhibits a similar trend to those in our galaxies—characterized by a rapid acceleration phase transitioning to an asymptotic velocity plateau.

To further assess this similarity, we fit each curve by varying v_∞ and β while fixing r₅₀ to the observed values. An illustrative example is shown for J1416+1223 (orange dashed line), where a β law with v_∞ = 1087 km s⁻¹, r₅₀ = 0.27 kpc, and β = 0.8 closely traces the inferred model curve (the top orange solid line). The right panel of Figure 5 presents the results for the full sample. We find that nearly all v_cl(r) profiles are well reproduced using β = 0.5–1.5, consistent with the resolved outflow structure observed in M82. These are different from the alternative forms of v_cl(r) proposed in other studies, e.g., M. Huberty et al. (2024) adopt . We highlight the need for future large samples of spatially resolved observations to definitively constrain the underlying outflow velocity structure.

In Figure 6, we present the radial distributions of η_M,cool(r) and η_M,hot(r), the outflow mass-loading factors for the cool and hot phases, separately. They are calculated from the radial-dependent mass-loading factor (Figure 3) divided by the SFR of the galaxy. The majority of the galaxies show a monotonically increasing η_M,cool(r) and a decreasing η_M,hot(r) with radius. There are a few exceptions that show completely the reverse trend (J1414+0540, J1448–0110, and J1521+0759 in both panels, and J0150+1308 and J0942+3547 in the right panel). We do not find special properties of them versus other galaxies, except they have relatively high fitted values for M_cl,0 (≳10⁵ M_⊙), but we note our adopted single value for M_cl,0 is likely not realistic (see the discussion in Section 5.2). It is likely that those outflows fail to be accelerated due to their low mass loading. This pattern has also been predicted in previous simulations (e.g., E. E. Schneider et al. 2020), which we discuss in Section 5. In both panels, we also see that more massive galaxies tend to have smaller loading factors, which is consistent with known scaling relationships of outflows (e.g., T. M. Heckman et al. 2015).

In Figure 7, we present the variation of the energy-loading factor over radius (η_E(r)), which is the energy outflow rate normalized by the the SF energy injection rate (). is computed using Starburst99 models (C. Leitherer et al. 1999), yielding = 4.3 × 10⁴¹ × SFR ergs s⁻¹. In the top two panels, we show η_E(r) for the cool and hot phases, separately. In all galaxies, we find η_E,cool(r) (η_E,hot(r)) increases (decreases) with radius. This is again expected given energy transfer from the hot to cool phases. In the bottom panel, we sum the energy-loading factors from the two phases to check the total energy of the system. We find the multiphase wind + outflow component loses substantial amounts (10%–90%) of its initial energy while it moves out of the galaxies. In Figure 8, we compare the energy loss in the two most likely scenarios, i.e., in TRML (L_TRML) or through radiative cooling in the volume-filling hot phase (L_hot). It is clear that the former dominates. This result suggests that turbulent mixing between the hot wind and entrained cool clouds provides the primary channel for radiative energy loss in multiphase galactic winds.

Zoom In Zoom Out Reset image size

Overall, the interactions between the cool and hot phases in galactic winds significantly affect the radial distributions of both components. While most of our results are consistent with empirical expectations, some rare cases show unanticipated results. We do not pursue these cases further to avoid overinterpretation and we will return to these questions in Paper II. In addition, the substantial mass, momentum, and energy exchanges between the hot wind to the cool clouds also underscores that galactic-wind models and simulations need to capture both phases. Otherwise, inferences from single-phase treatments can be biased.

Galactic winds are of great interest because they are a primary mechanism driving feedback that shapes galaxy evolution. Over the last decade, significant progress has been made in understanding both the observations and the multiphase physics of these outflows. These include multiwavelength observational studies of galactic winds, as well as analytic and hydrodynamic models spanning a wide range of resolutions and domain sizes (see the reviews by S. Veilleux et al. 2020; T. A. Thompson & T. M. Heckman 2024).

Specifically, the prototypical starburst galaxy M82 (SFR = 8 M_⊙ yr⁻¹, M_⋆ = 10¹⁰ M_⊙, and r₅₀ = 300 pc) has accumulated extensive studies, owing to its proximity and well-defined biconical outflows (e.g., P. L. Shopbell & J. Bland-Hawthorn 1998; D. K. Strickland & T. M. Heckman 2009; A. K. Leroy et al. 2015; P. Martini et al. 2018). In observations, we have measured the radial distribution of cool outflow cloud properties using Subaru/FOCAS (X. Xu et al. 2023a), and overlaid our results as black lines with gray bands (1σ error) in Figures 5–7. First of all, the observed v_cl in M82 follows our model trends and closely tracks one of our galaxies, J1112+5503, which is also a massive, highly star-forming system (SFR ∼ 40 M_⊙ yr⁻¹ and M_⋆ = 10^9.6 M_⊙). For the mass- and energy-loading factors (Figures 6 and 7), the M82 measurements are consistent with either flat radial evolution or with declines, while our models show similar behavior for massive galaxies.

Furthermore, S. Lopez et al. (2025) study the properties of cool outflow clouds in M82 from Hα-emitting structures using HST/Advanced Camera for Surveys imaging from M. Mutchler et al. (2007). They report that individual cool clouds commonly exhibit elongated elliptical or arc-like morphologies rather than the cometary structures predicted in some simulations (e.g., M. Gronke & S. P. Oh 2018; E. E. Schneider et al. 2020; FB22). Their inferred cloud column density (N_H,cl) is ∼10²⁰–10²¹ cm⁻², matching well with our galaxies (column (7) in Table 1). Their cloud masses, M_cl, lie in the range 10⁴–10⁵ M_⊙, consistent with the lower bounds in our best-fit models (column (10)). Their resolved clouds reside at radii of ∼0.5–2 kpc (∼2–7 × r₅₀).

Observational studies of the hot and very hot phases of galactic winds are much rarer, typically detected via extended, diffuse X-ray emission (e.g., D. K. Strickland & T. M. Heckman 2009; J. Liu et al. 2012; D. Zhang et al. 2014; L. A. Lopez et al. 2020; K. Iwasawa et al. 2023). In M82’s central starburst (r ≲ r₅₀), hard X-ray spectra are dominated by He- and H-like Fe at T ∼ 10⁸ K. The inferred terminal velocity is 1400–2200 km s⁻¹ and the hot-phase mass loading is ∼0.2–0.6 (e.g., D. K. Strickland & T. M. Heckman 2009). These values again are consistent with our modeled results for massive galaxies (e.g., right panel of Figure 6). Using soft X-ray data, E. Boettcher & E. Hodges-Kluck (2024) analyze deep archival Reflection Grating Spectrometer observations from XMM-Newton and detect clear emission from O viii (0.65 and 0.77 keV), Ne x (1.02 keV), and Mg xii (1.47 keV). They infer a hot-wind phase with T ∼ 3 × 10⁶ K and velocity ≳ 2000 km s⁻¹ (from O viii), consistent with the hard X-ray estimates and with the model predictions for massive galaxies in our sample.

Several simulations directly model properties of the hot wind. E. E. Schneider et al. (2020) present high-resolution (<5 pc) simulations of a multiphase outflow from a disk galaxy with a starburst radius of 1000 pc, M_⋆ = 10¹⁰ M_⊙, and SFR = 5–20 M_⊙ yr⁻¹. They highlight that mixing between hot (in their case, T > 5 × 10⁵ K) and cool (T < 2 × 10⁴ K) gas efficiently transfers momentum at all radii. In their high-SFR state (SFR = 20 M_⊙ yr⁻¹) at 35 Myr, η_M,hot increases from 0.04 to 0.06 between 2 and 8 r₅₀, while η_M,cool decreases from 0.04 to 0.01 over the same range. These values are lower than the ones for most of our galaxies, but they approach the lower bounds in Figure 6, where we indeed find cases with decreasing η_M,cool, increasing η_M,hot, and both with values < 0.1. This likely reflects a regime in which the hot wind fails to launch strong cool outflows.

By contrast, other findings predict consistent hot-phase properties with most of our galaxies. For example, X.-F. Li et al. (2025) simulate a series of starbursts with SFR = 1.5–30 M_⊙ yr⁻¹, M_⋆ = 1.6–3.3 × 10⁸ M_⊙, and the majority of stars forming within the central 1000 pc. They find η_M,hot ∼ 0.2–0.4 at 10 r₅₀ and hot-wind velocities of 800–1400 km s⁻¹. In addition, D. D. Nguyen & T. A. Thompson (2021) develop an analytic hot-wind model that includes mass loading and nonspherical divergence, calibrate it to simulations, and apply it to M82. They assume SFR = 10 M_⊙ yr⁻¹ and r₅₀ = 300 pc. They find a total mass-loading factor of about one and predict an asymptotic hot-wind velocity of ∼1000 km s⁻¹. They also emphasize that geometry and mass loading can commonly lower v_wind and η_M,hot relative to non-mass-loaded spherical expectations.

Overall, our work combines observations and models of the cool and hot phases across a wide range of galaxy properties, and the resulting trends agree with prior multiphase, multiscale studies of M82.

Our current methodology and the model have several caveats and some of them are tightly related to our future directions. First of all, the model assumes a single initial cloud mass, M_cl,0, for all cool outflow clouds; this simplification precludes recovering a realistic cloud-mass distribution. Then the energy-loading factors (η_E)—particularly for the hot phase—are not included as free parameters in the present implementation. In CC85 framework, hot-phase properties scale with the ratio η_E/η_M, so allowing η_E to vary is essential for completeness. In the follow-up Paper II, we will modify the methodology and model hierarchy to address these two points, which will help us understand why several galaxies are not well fit under the current settings.

Another longstanding question is what is the best way to extract the hydrogen column density profile (i.e., dN_H/dv) from rest-frame UV absorption lines. We adopt the FB22 model to fit dN_H/dv derived using the PCM from CLASSY data (X. Xu et al. 2022). More recent studies have shown that the PCM likely yields underestimated column densities. For example, M. Huberty et al. (2024) and C. A. Carr et al. (2025) suggest that the underestimation is due to the finite instrumental resolution and the resulting smoothed spectral data. On the contrary, R. M. Jennings et al. (2025) suggest that such underestimation is potentially caused by missed high-column-density sight lines that are optically thick, which are invisible in UV spectra. While the exact biases introduced by the PCM—if any—remain under active investigation, the methods presented in this paper are applicable to any dN_H/dv distribution, provided the outflow velocity field is monotonic.

Finally, the current FB22 model does not account for the CGM component. As the galactic wind expands, it can sweep up substantial amounts of ambient CGM material, potentially altering the mass-loading factors. Future comparisons with the results from cosmological simulations, such as FIRE-2 (P. F. Hopkins et al. 2018), will help assess the magnitude of this effect.