CLASSY III: The Properties of Starburst-Driven Warm Ionized Outflows

We report the results of analyses of galactic outflows in a sample of 45 low-redshift starburst galaxies in the COS Legacy Archive Spectroscopic SurveY (CLASSY), augmented by five additional similar starbursts with COS data. The outflows are traced by blueshifted absorption-lines of metals spanning a wide range of ionization potential. The high quality and broad spectral coverage of CLASSY data enable us to disentangle the absorption due to the static ISM from that due to outflows. We further use different line multiplets and doublets to determine the covering fraction, column density, and ionization state as a function of velocity for each outflow. We measure the outflow's mean velocity and velocity width, and find that both correlate in a highly significant way with the star-formation rate, galaxy mass, and circular velocity over ranges of four orders-of-magnitude for the first two properties. We also estimate outflow rates of metals, mass, momentum, and kinetic energy. We find that, at most, only about 20% of silicon created and ejected by supernovae in the starburst is carried in the warm phase we observe. The outflows' mass-loading factor increases steeply and inversely with both circular and outflow velocity (log-log slope $\sim$ -1.6), and reaches $\sim 10$ for dwarf galaxies. We find that the outflows typically carry about 10 to 100% of the momentum injected by massive stars and about 1 to 20% of the kinetic energy. We show that these results place interesting constraints on, and new insights into, models and simulations of galactic winds.


INTRODUCTION
We live in a time of challenges and opportunities in the quest to understand the evolution of galaxies. We have a very successful theory for the development of the large-scale structure of the dark matter scaffolding in which the galaxies form and grow (e.g., Wechsler & Tinker 2018). We also know the overall cosmic history of the rate of galaxy buildup through measurements of the star-formation rate (SFR) per unit co-moving volume element (e.g., Madau & Dickinson 2014) and through measurements of the evolution of the cosmic inventory of baryons (e.g., Péroux & Howk 2020).
In the current paradigm of galaxy evolution (e.g., Somerville & Davé 2015;Naab & Ostriker 2017), baryons flow with the dark matter on large scales, and some are incorporated into the halos. Unlike the dark matter, the baryons can lose energy by radiation, and sink deeper into the potential well defined by the dark matter. In the simplest picture, this inflow is halted by centrifugal forces reflecting conservation of angular momentum. Stars form within the central regions of these disks. Galaxies continue to grow over billions of years, primarily through continuing accretion of gas from the cosmic web, and secondarily through mergers with other dark matter halos and their baryonic contents.
While this picture is simple and compelling, it does not account quantitatively for even the most basic properties of the baryonic content of galaxies. Some of the key unsolved problems are (e.g., Somerville & Davé 2015;Behroozi et al. 2019;Maiolino & Mannucci 2019;McGaugh et al. 2000): 1. Why are only about 10% of the baryons accreted (as gas) into the dark matter halos incorporated into stars, and why is this efficiency largely independent of redshift (z)?
2. Why does this efficiency reach its peak value over a relatively narrow range in dark matter halo masses of ∼ 10 12 M , and why is this value roughly constant with z?
3. Why is there such a tight correlation between a galaxy's stellar mass (M ) and its SFR at a given epoch (the "star-forming main sequence")?
4. Why is there such small scatter in the correlation between M and the galaxy's chemical composition (metallicity) at a given epoch? 5. Why is the scatter so small in the relationships between M , internal velocity dispersion and/or rotation speed, and radius in galaxies at a given z, and why do they evolve with z?
6. How does the intergalactic medium (IGM) get enriched with metals?
In all current theoretical models and numerical simulations of galaxies, these questions are dealt with through the rubric of "feedback": the effects of the return of energy, momentum, and heavy elements from massive stars and black holes on the surrounding gas (e.g., Somerville & Davé 2015;Naab & Ostriker 2017). The tightness of the scaling relations for galaxies (questions 3, 4, and 5 above), requires true two-way feedback that leads to self-regulating processes.
As noted above, feedback can be provided by either stars or supermassive black holes. In this paper, we focus on the former. Stellar feedback is dominated by massive stars, and is supplied in the form of radiation and stellar ejecta. For a young stellar population, the kinetic energy and momentum are primarily supplied by a combination of stellar winds from hot, massive stars and through core-collapse supernovae (e.g., Leitherer et al. 1999). In order to affect the structure and content of a galaxy, either the momentum or kinetic energy provided by these stars must couple to the gas supply of the galaxy.
Undoubtedly, the most spectacular manifestations of feedback from populations of massive stars are global-scale galactic winds (e.g., Heckman & Thompson 2017;Veilleux et al. 2020). These winds play a crucial role in the evolution of galaxies and the IGM (e.g., Somerville & Davé 2015;Naab & Ostriker 2017). In particular, the selective loss of gas and metals from the shallower potential wells of low-mass dark matter halos is believed to be responsible for both shaping the low-mass end of the galaxy stellar mass function and for establishing the mass-metallicity relation (questions 2, and 4 above). By carrying away low-angular-momentum gas, they also shaped the mass-radius relation (question 5). These outflows heated and polluted the circum-galactic medium (CGM) and IGM with metals and may have suppressed the accretion of gas passing from the CGM into the star-forming disk (questions 1 and 6).
Simply put, the evolution of galaxies cannot be understood without first understanding galactic winds. Given their importance, and given both the large amounts of data that have been collected and the increasing quality of numerical simulations, it is perhaps surprising that we still have a very incomplete understanding of the processes that create outflowing gas, and of the impact the outflow has on the galaxy that launches it.
It is crucial to emphasize that feedback processes associated with galactic winds cannot be spatially-resolved in numerical cosmological simulations. This problem has been long-recognized (e.g., White & Frenk 1991;Katz et al. 1996;Springel & Hernquist 2003;Hopkins et al. 2014). Instead the feedback processes are implemented numerically using "subgrid physics" (recipes, see e.g., Somerville & Davé 2015;Naab & Ostriker 2017). These recipes often depend on things like the SFR and mass of the galaxy. It is clear that observations of feedback in-action are essential in guiding these choices and revealing the actual dependences on galaxy properties. Such data are also required to test models which attempt to simulate galactic winds with high enough spatial resolution to allow more ab initio calculation of the relevant physics (e.g., Schneider et al. 2018;. To date, the bulk of the data on winds across cosmic time have come from analysis of interstellar absorption-lines that trace outflowing cool or warm gas through the blue-shifted absorption-lines it produces (Heckman et al. 2000;Shapley et al. 2003;Rupke et al. 2005;Martin 2005; Grimes et al. 2009;Sato et al. 2009;Weiner et al. 2009;Rubin et al. 2010;Steidel et al. 2010;Chen et al. 2010;Erb et al. 2012;Kornei et al. 2012;Martin et al. 2012;Bordoloi et al. 2014;Rubin et al. 2014;Heckman et al. 2015;Zhu et al. 2015;Chisholm et al. 2015;Heckman & Borthakur 2016;Chisholm et al. 2016a;Sugahara et al. 2017;Steidel et al. 2018;Chisholm et al. 2018;Sugahara et al. 2019).
Clearly, there have been many prior investigations. In this paper, we seek to significantly improve the usefulness of such data in two respects. First, we use the COS Legacy Archive Spectroscopy SurveY (CLASSY) atlas ; hereafter, Paper I), which is a data set specifically designed to span a vast range in the most fundamental galaxy properties: stellar mass (and hence galaxy circular velocity), star-formation rate, and metallicity (see details in Section 2). Reaching low mass is especially important since feedback effects should be stronger in these shallow potential wells. This makes CLASSY ideal for testing how the fundamental outflow properties (outflow velocities, column densities, ionization state, and metal, mass, momentum, and kinetic energy outflow rates) depend on the galaxy properties, thereby providing a crucial test of the theoretical models and simulations. Second, the CLASSY data, by design, are high signalto-noise ratio (SNR) spectra and cover a wide wavelength range in the UV that encompasses many interstellar lines that span wide ranges in ionization state and optical depth. Both points enable us to analyze the data more rigorously than be-fore, leading to more robust measurements of outflow properties.
The structure of the paper is as follows. In Section 2, we introduce the CLASSY project and briefly describe the observations. In Section 3, we go through various data reduction processes. In Section 4, we present analyses to isolate the blue-shifted absorption lines for galactic outflows, and we also discuss the ancillary parameters. In Section 5, we present the major results for the observed outflows, including covering fraction, column density, and ionization state as a function of velocity, and we also derive outflow rates of metals, mass, momentum, and kinetic energy. Furthermore, we compare the derived outflow properties to various host galaxy characteristics. Finally, in Section 6, we discuss and compare our results to common models for galactic winds, as well as semi-analytic models and numerical simulations. We summarize the paper in Section 7.

OBSERVATIONS
CLASSY is a Hubble Space Telescope (HST) treasury program (GO: 15840, PI: Berg), which provides the first highresolution high-SNR restframe far-ultraviolet (FUV) spectral catalog of 45 local star-forming galaxies (0.002 < z < 0.182). These galaxies were selected to span a wide range of important galaxy properties, including stellar mass (log M ∼ 6 − 10 M ), SFR (∼ 0.01 − 100 M yr −1 ), metallicity (12+log(O/H)∼ 7 − 9), and electron density (n e ∼ 10 1 − 10 3 cm −3 ). For each galaxy, CLASSY completes its FUV wavelength coverage (1200 -2000 Å observed frame) by utilizing the G130M+G160M+G185M/G225M gratings of HST/COS. Overall, CLASSY combines 135 orbits of new HST data with 177 orbits of archival HST data to complete the first atlas of high-quality restframe FUV spectra of the proposed 45 galaxies. We define the CLASSY data hereafter as this combined dataset. We refer readers to Paper I for the detailed sample selection, observations, and basic properties of these galaxies.

DATA REDUCTION
After the observations, all data were reduced locally using the COS data-reduction package CalCOS v.3.3.10 1 . These include both new data from CLASSY itself (GO: 15840) and all archival data that were included in CLASSY. Therefore, the whole CLASSY dataset was reduced and processed in a self-consistent way. The details of the data reduction have been presented in Paper I and Paper II, including spectra extraction, wavelength calibration, and vignetting.
Given the final reduced and co-added spectra for each galaxy (Paper I), we analyze the galactic outflow properties from various absorption and associated emission lines in this paper. Several additional steps in the reduction of the data are necessary, and are discussed in this section. In Section 3.1, we discuss the fits of the stellar continuum for each galaxy, which are used to remove the starlight contamination of the spectra. In Section 3.2, we discuss other systematic effects in the analyses of outflows.
To enlarge the sample size at the highest star-formation rates, we have added five Lyman Break Analog (LBA) galaxies from Heckman et al. (2015) that were not already in the CLASSY sample but also had HST/COS observations. We have checked that these LBAs satisfy all selection criteria of the CLASSY sample. We processed and analyzed these data in exactly the same way as the CLASSY data. This gives us a total sample size of 50 galaxies. HST/COS G130M and G160M gratings have the original spectral resolutions R ∼ 20,000, which is more than necessary for our outflow analyses. Therefore, for all spectra, we re-sample them into bins of 0.18 Å (R ∼ 6000 -10000 from the blue to red end) to gain higher S/N.

Starlight Normalization
By assuming that the observed spectra are combinations of multiple bursts of single-age, single-metallicity stellar populations, one can fit the stellar continuum of galaxies by linear combinations of stellar models, e.g., from Starburst99 (Leitherer et al. 1999). We do so by following the same methodology laid out in Chisholm et al. (2019). Then, for each galaxy, we normalize the spectra by the best-fit stellar continuum.

Other Systematic Effects
To get robust measurements of the outflow properties from the absorption lines, we need to take account of multiple systematic effects or contamination: 1) the static ISM component that is centered at v = 0 km s −1 . This component represents the ISM gas that is not accelerated by the galactic winds, but it can blend with the lower velocity portions of the outflows; 2) the HST/COS line-spread-functions (LSFs), which describe the light distribution at the focal plane as a function of wavelength in response to the light source. This can slightly broaden and reshape the observed absorption troughs; and 3) the infilling effects from corresponding resonantly-scattered emission lines. We quantify each of these points in our analyses in Section 4 below.

BASIC ANALYSES
We begin with presenting a brief justification of the major methods adopted in Section 4.1. Then in Section 4.2, we show the double-Gaussian fits to the observed outflow absorption troughs. In Section 4.3, we discuss the infilling of Figure 1. Comparisons of the line spread functions (LSFs) for galaxy J1148+2546. The point-source LSF is shown in red, which is for the HST/COS G130M grating given the center wavelength of 1222 Å and at Lifetime Position 4 (LP4). The Gaussian representing the NUV light profile in the dispersion direction of this galaxy is shown in blue. The approximate LSF used in our fits is shown in black, which is the convolution between the blue and red curves. This galaxy has a relatively large NUV size in our sample, so the differences between the two LSFs are noticeable. See discussion in Section 4.2.
absorption troughs by corresponding emission lines. Finally, in Section 4.4, we discuss various ancillary parameters that are adopted in the rest of the paper.

Justifications of Methodology
For each galaxy in our sample, the final reduced, co-added, and starlight subtracted spectra cover ∼ 1200 Å -2000 Å in the observed frame. In this region, various lines from galactic outflows are detected as absorption troughs, including from low-ionization transitions, e.g., O I λ1302, C II λ1334, and Si II multiplet (λ1190, 1193, 1260, 1304, and 1526), and from higher-ionization transitions, e.g., Si III λ1206, and Si IV λλ1393, 1402. We focus on these lines in our fitting method described in this section. In Table 1, we list important atomic information for these lines.
To determine the basic properties of the outflows, we first fit the observed absorption troughs (Section 4.2). In the meantime, we need to take account of the spectral linespread function (LSF) which is convolved with the intrinsic absorption-line profile to produce the observed profile. The LSF has contributions from both the COS optics and the spatial distribution of the UV continuum in the COS aperture. Our approach is to fit the observed profiles by using a simple analytic form for the intrinsic line profile, which has been convolved with our calculated LSF for each galaxy. We will describe this in more detail below. Here we note that we have taken a Gaussian to describe the intrinsic absorptionline profiles for both a component associated with the static ISM and one associated with the outflow. The choice of a Gaussian is motivated by several considerations. First, it is a simple analytic function, unlike the Voigt profile (Draine Figure 2. Examples of double-Gaussian fits and F-tests to the absorption troughs. The normalized fluxes are shown in black histograms. For double-Gaussian fits, the profiles for the outflow and static ISM components are shown in blue and green dashed lines, respectively, and their summation is shown in red. For single-Gaussian fits, the models are shown in purple solid lines. All models have already considered the LSF effects discussed in Section 4.2. The χ 2 values for single-and double-Gaussian fits are listed at the bottom left corner of each panel in purple and red, respectively. The fitting ranges are within the two vertical gray dashed lines. Left: A case when the absorption trough passed the F-test. This Si II λ1260 trough has a significant blue shift. Thus, the second Gaussian (in blue) is necessary to fit the trough and represent the outflow. Right: A case when the trough failed the F-test. In this case, we don't gain significant improvement by introducing the second Gaussian, so we label this trough as "no outflow". See details in Section 4.2.   (Kramida et al. 2018).
(5). Energies from lower to higher levels in units of eV. 2011). Second, it provides an excellent fit to the data (as we will show).

Double-Gaussian Fits of the Absorption Troughs
The theoretical LSF provided in HST/COS website 2 is good for a point source, but CLASSY galaxies are largely 2 https://www.stsci.edu/hst/instrumentation/cos/performance/ spectral-resolution resolved. Therefore, we need to construct non-point source LSF for each galaxy, separately. The steps along with the double-Gaussian fitting process are as follows.
1. First of all, we would like to generate a LSF for each galaxy based on the COS LSF for a point source (LSF 0 ) and the spatial distribution of the FUV continuum in the COS aperture in the dispersion direction. For the latter, we consider the galaxy size that is measured from HST/COS NUV acquisition images (FWHM uv , see Table 2 in Paper I). An example is shown in Figure 1). We assume the galaxy has Gaussian profile with FWHM uv (G uv , in blue), which is sufficient for our analyses. Then we convolve LSF 0 (in red) with G uv , which results in an approximate nonpoint source LSF (LSF uv , in black) for a given galaxy. This LSF uv is then used in the double-Gaussian fittings below to properly account for the LSF of each galaxy.
2. We then use a double-Gaussians model to fit the observed absorption troughs adopting the fitting routine mpfit (Markwardt 2009). Two examples are shown in Figure 2. Instead of using the standard Gaussian profile, we convolve it with LSF uv measured in step 1 to take into account the effects of LSFs (hereafter, we use G * 1 and G * 2 for the two convolved Gaussian profiles). G * 1 has a fixed velocity center at v = 0 km s −1 , which accounts for the static ISM component (green lines in Figure 2). G * 2 has a velocity center < 0 km s −1 that represents the outflow component (blue lines in Figure  2).
3. For each fitted absorption trough, to check if G * 2 is necessary (i.e., if there exists an outflow component), we conduct an F-test: where χ 2 1 and χ 2 2 are the chi-squares from a single-Gaussian model (fits of the trough by only G * 1 ) and a double-Gaussian model (fits of the trough by G * 1 and G * 2 ), respectively. p 1 and p 2 are the number of free parameters in single-and double-Gaussian models, respectively, and n is the total number of bins of the fitted trough. We then compare this calculated F value with the theoretical one from F-distribution table given significance level α = 0.05 and degree of freedom of (p 2p 1 , np 2 ). If the fitted F value is greater than the theoretical one, we reject the null hypothesis (i.e., model 2 does not provide a significant better fit than model 1). This indicates that the inclusion of G * 2 is necessary to fit the observed trough. Therefore, we treat this trough hosting blueshifted outflow(s). In Figure 2, we show examples of passing/failing the F-test in the left/right panels, respectively.
For most galaxies, we find that the effects of the LSFs on the outflow components are relatively small. This is because 1) the LSFs do not alter the measured velocity centers of the fitted outflow components, and 2) the FWHMs of LSF uv are usually 100 km s −1 , but the FWHMs of the outflow components are usually > 250 km s −1 . On the contrary, for the narrow static ISM component and galaxies with narrow outflow troughs (< 100 km s −1 ), the HST/COS LSFs have significantly broadened their FWHMs. Therefore, in these cases, our fitting method discussed above is necessary for quantifying the outflow's FWHM (FWHM out ): where FWHM all is the fitted FWHM for a certain trough and FWHM LSF is the FWHM of LSF uv for that galaxy.
For each galaxy, we adopt this method to fit the absorption lines from O I λ1302, C II λ1334, Si II multiplet (λ1190, 1193, 1260, 1304, and 1526), Si III λ1206, and Si IV λλ1393, 1402, separately. If one individual line falls into a chip gap or is contaminated by Galactic lines (e.g., Si III can be affected by Galactic Lyα), we exclude them from the fitting. For each galaxy, if more than half of the fitted troughs for the different absorption-lines pass the F-test, we label this galaxy as "hosting outflows". Otherwise, if one galaxy has less than half of its troughs that passes the F-test, "no outflow" is labelled.
One galaxy's outflow velocity is then calculated from the median value of central velocities (of fitted G * 2 ) from all troughs that have passed the F-test. Similarly, the galaxy's outflow full-width-half-maximum (FWHM out ) is derived from the median value of FWHM from all troughs that pass the F-test. The corresponding errors of V out (or FWHM out ) are estimated from the standard deviations of V out (or FWHM out ) from all passed lines. Note that for these "no outflow" galaxies, they could still host very low velocity outflows (V out FWHM ISM ). However, since we cannot disentangle them from the static ISM component, their V out and FWHM out are not measurable.
Besides the median values, we also have examined the consistency of the derived outflow velocity and FWHM among the different transitions. Two examples are shown in Figure  4. We find that the values of outflow velocities are quite consistent, but there is significant scatter in the values of FWHM. The panel showing FWHM for the two Si II transitions is particularly instructive. It shows a systematic trend for the FWHM from Si II λ1304 (which is the most optically thin Si II line we measure) to be narrower than that of Si II λ1260 (which is the most optically thick Si II line that we measure). We will discuss the implications of this in Sections 5.1 below.
Overall, 43 out of 50 galaxies (∼ 86%) in our combined sample are labeled as "hosting outflows". This indicates that galactic outflows are common in the low-redshift starburst galaxies in the CLASSY sample. The distribution of V out and FWHM out are shown in Figure 3, and the values are presented Table 5.

Effects of Infilling
As discussed in previous publications (e.g., Prochaska et al. 2011;Scarlata & Panagia 2015;Alexandroff et al. 2015;Carr et al. 2018;Wang et al. 2020;Carr et al. 2021), the absorption of a photon from the ground state may be followed by emission (resonance scattering). While only the gas directly along the line-of-sight to the UV continuum will produce absorption, the emission will come from all the gas within the size of COS aperture. This line emission will infill the absorption trough and could significantly affect our fitted models. Therefore, it is necessary to check the effects of infilling during our fitting processes.
The infilling is difficult to disentangle given only the absorption lines themselves, but fortunately, many of these lines have associated fluorescent emission lines. These are the fine structure transitions that share the same upper energy level as the associated resonance lines, e.g., Si II* λ1265 for Si II λ1260. Since the fluorescent lines are well displaced in wavelength from the resonance lines, their properties can be easily determined. Furthermore, the physical properties of the infilled emission line is related to that of fluorescent emission lines in various ways. For example, the infilling emission of the Si II λ1260 would cover the same velocity range as Si II* λ1265, and the line strength ratio of 1260/1265 depends on: 1) their Einstein-A coefficient ratios; and 2) how many times the photon has been scattered (represented by the optical depth of the trough, see, e.g., Scarlata & Panagia 2015). We mainly test infilling effects on the Si II multiplet since they are the major lines affected and commonly used in our analyses later (see Section 5). The other major lines are Si IV λλ1393, 1402, but there are no associated fluorescent transitions for them. The steps are as follows: 1. For each Si II resonance line, we first fit the corresponding fluorescent emission line adopting one convolved Gaussian (G * , see solid orange lines in Figure  5).
2. Then we choose the ten objects that have the largest |EW(Si II*)/EW(Si II)|, since we expect that these galaxies are affected the most by infilling. We conduct a triple-Gaussian fit to their Si II absorption troughs. Besides the double Gaussians (in absorption) discussed in Section 4.2, we add a third Gaussian (G * 3 , in emission) to represent the infilling component. We fix the velocity center and width of G * 3 as the same as its associated fluorescent line (fitted in step 1) but allow the amplitude of G * 3 to vary between 0 and a constant. For Si II λ1260, this constant equals A 1260 /A 1265 × H 1260 0.85H 1260 (e.g., Scarlata & Panagia 2015), where A is the Einstein A coefficient and H is the amplitude of the emission line.
3. Finally, we compare the results from the double-Gaussian fits and triple-Gaussian fits by measuring three parameters, i.e., velocity center, FWHM, and minimum flux of their outflow trough (G * 2 ).
Overall, for these ten galaxies, these measured parameters for G * 2 only have minimal differences (< 5%) between the double-Gaussian and triple-Gaussian fits. Thus, we conclude that in our sample, the infilling is negligible for the outflow absorption troughs from Si II multiplet (e.g., Alexandroff et al. 2015). Nonetheless, as shown in the right panel of Figure 5, infilling can more significantly affect the properties of the fit to the static ISM component (centered at v = 0 km s −1 ), which would alter the derived covering fraction, column density, etc., for the ISM component. These conclusions are consistent with extensive studies presented in previous publications of metal absorption lines in galaxies (e.g., Prochaska et al. 2011;Martin et al. 2012;Erb et al. 2012;Alexandroff et al. 2015).

Ancillary Parameters
Various galaxy properties are used in the rest of the paper. The measurements of them are discussed in Section 4 of Paper I and listed in their Tables 5 and 6. We take these values and propagate their errors when they are used in our analyses. These include: 1) gas-phase metallicity (12+log(O/H)) and the errors, which are derived using the direct method; 2) M , SFR, and dust extinction (E(B-V)), and the corresponding errors, which are derived from the SED fitting to the GALEX+SDSS data; 3) redshift of galaxies, which are derived from fitting the optical emission lines (mainly from SDSS spectra, Mingozzi et al. in prep.).
We also calculate two other ancillary parameters needed for our analyses, including the galaxy circular velocity (V cir ) and r 50 . Following Heckman et al. (2015), we calculate V cir based on the galaxy stellar mass. This adopts the empirical calibration in Simons et al. (2015) which was derived using spatially-resolved maps of the kinematics of the gas in low-redshift emission-line galaxies. Following Simons et al. (2015), we define V cir = √ 2 (V 2 rot + 2σ 2 ) 1/2 , where V rot and σ are the rotation speed and mean value of the line-of-sight velocity dispersion, respectively. Finally, we adopt the bestfit relation from Simons et al. (2015) as log V cir = 0.29 M -0.78, with root-mean-square (RMS) residuals of 0.1 dex RMS in V cir .
For each galaxy, we measure the half-light radius (r 50 ) from its HST/COS NUV acquisition image. We accumulate the net photon counts within a certain radius from the galaxy center until it reaches half of the total source counts. However, some galaxies are more spatially extended than the unvignetted region of the acquisition images, and their cumulative counts vs. radius distributions are not flat at large radii. Therefore, for the scaling relations that we discuss in Section 5.4, we adopt the r 50 values from COS for galaxies with r 50 smaller than 0.4 . These galaxies are compact enough, so the vignetting effects are small. For other galaxies that have COS r 50 > 0.4 , we instead take their SDSS u-band sizes from the SDSS archive, which conducted exponential fits convolved with the measured PSF to the SDSS images.
For three objects (J1129+2034, 1314+3452, J1444+4237) in our sample, the SDSS u-band sizes from the archive only refer to a bright star-formation knot instead of the entire galaxy. Therefore, we remeasure their sizes using the SDSS u-band image and performing the method of aperture photometry discussed above. Furthermore, there are five galaxies (J0036-3333, J0127-0619, J0144+0453, J0405-3648, J0337-0502) that do not have SDSS u-band images. We instead measure sizes from their optical images that are close to the u-band, including Dark Energy Survey (DES) g-band and Pan-STARRS g-band images.
For galaxies with r 50 > 1.5 , the COS aperture does not cover the majority of massive stellar population in the galaxy. This means that their observed absorption-line outflows may be less related to the global properties of the galaxy but in- stead more related to bright local properties (e.g. individual massive young clusters). Therefore, we will label these galaxies differently (in gray triangles) in all correlation figures in Section 5. For each galaxy, the final adopted V cir and r 50 values along with other important ancillary parameters used in this paper are listed in Table 6. More detailed studies of the aperture effects of CLASSY galaxy properties are presented in Arellano-Córdova et al. submitted.

RESULTS
In this section, we present detailed analyses and the corresponding results for galaxies in the CLASSY sample. We start with getting robust estimates of the column density and covering fraction (C f ) of outflows from the observed Si II and Si IV absorption lines in Section 5.1. Then we study the dust depletion of metals in these galaxies from stacked spectra in Section 5.2. After that, we estimate the total hydrogen column density (N H ) of outflows from CLOUDY models (Ferland et al. 2017) in Section 5.3. Finally, we present various scaling relationships related to outflow kinematics, feedback effects, and galaxy properties in Section 5.4.

Column Density and Covering Fraction of Si II and Si IV Lines
Galactic outflows have been found to only partially cover the ionizing source (e.g., Heckman et al. 2000;Rupke et al. 2005;Martin & Bouché 2009;Chisholm et al. 2016a). In the case of covering fraction (C f ) < 1.0, the column density (N ion ) derived from the absorption lines assuming an apparent optical depth (AOD) can only be viewed as a lower limit. Therefore, to get robust N ion and C f measurements, we adopt the partial-covering (PC) models that were commonly used in analyzing both quasar outflows (e.g., Arav et al. 2005;Xu et al. 2018) and galactic outflows (e.g., Rupke et al. 2005;Martin & Bouché 2009;Chisholm et al. 2016a;Rivera-Thorsen et al. 2015;. These models assume that only a portion of the UV continuum source (in our case, starburst region) is covered by foreground absorbing/outflowing gas. C f and the optical depth (τ ) are usually degenerate since they can both affect the depths of a absorption line. However, the degeneracy can be broken given the useful information from doublet and multiplet transitions as follows (e.g., Arav et al. 2005).
For absorption lines from a doublet transition (e.g., Si IV λλ1393, 1402), we have: where I R (v) and I B (v) are the red and blue doublets' flux at velocity v, normalized by the stellar continuum (Section 3.1), respectively, and τ (v) is the optical depth of the red line (Si IV λ1402) that is velocity dependent. The weight w in the bottom equation equals f B λ B / f R λ R , where f is the oscillator strength and λ is the wavelength of the line (see Table 1). For common doublet lines, such as N V λλ1238, 1242, C IV λλ1548, 1550, and Si IV λλ1393, 1402, they have w = 2. For absorption lines from a multiplet transition (e.g., Si II multiplet), we similarly have a set of equations: where for Si II multiplet, k stands for Si II λ1190, 1193, 1260, 1304, and 1526. We define w k = f k λ k / f 1304 λ 1304 , and therefore w k = 2.7, 5.7, 12.7, 1.0, 1.7 for the five Si II lines, respectively. If one or more troughs from the multiplet is blended with other lines (e.g., Milky Way absorption lines) or in detector gaps, we exclude them from the equation set.
To solve Equations (3) or (4), we require the detections of absorption troughs from isolated doublet or multiplet transitions, respectively. Given the wavelength coverage of Figure 6. An example of the solutions of partial-covering equations and total hydrogen column densities (NH) for galaxy J0055-0021. Top-Left: The fitted Gaussian profiles for the four Si II outflow troughs are shown as the solid lines (i.e., G * 2 , Section 4.2). A total of four Si II troughs are clean (shown as different colors), while J0055-0021's Si II λ1526 trough is in a detector gap. Then we fit partial-covering models [see equation (4)] to these Gaussian profiles and the results are shown in the dashed lines. The errors of the Gaussian profiles are shown as dotted lines at the bottom of the troughs, and the fitting range is within the two gray lines. Note that the Si II λ1304 line is significantly narrower than the other three lines Top-Right: Best fitting covering fraction (Cf) for Si II from solving Equation (4). Bottom-Left: Best fitting optical depth over velocity (τ (v)) in logarithm for each Si II lines. We only show Si II λ1260's error bars to avoid crowding. Bottom-Right: The resulting average column density in the aperture per velocity for Si II, which is proportional to Cf(v)×τ (v) [Equation (4)]. In the latter 3 panels, the corresponding errors derived from mpfit are shown as vertical lines. This shows that the decline in column density away from the line center is due to drops in both optical depth and covering fraction. See details of the fitting process in Section 5.1. CLASSY data, the major lines discussed in the reminder of this paper are Si IV doublet and Si II multiplet. We didn't solve the absorption troughs from C IV λλ 1548, 1550 doublet. This is because the two troughs blend with each other, and the solutions are usually uncertain (e.g., Du et al. 2016). Furthermore, the C IV line is more likely to be affected by nebular emission and stellar P-Cygni absorption.
For each galaxy, our main goal is to solve for the N(Si IV) and N(Si II) in the galactic outflows. Therefore, we take fitted Gaussian models (i.e., G * 2 , see Section 4.2) which represents the outflow component. Note that, due to possible contamination from infilling at the systemic velocity (Section 4.3), it is difficult to get robust values for N(Si IV) and N(Si II) for the static ISM component (i.e., G * 1 ). We then adopt mpfit (Markwardt 2009) to solve τ (v) and C f (v) from Equations (3) and (4) for Si IV and Si II, respec-tively. An example for the fitted Si II troughs and the resulting C f (v), τ (v), and N ion are shown in Figure 6. For each galaxy, we also list the outflow mean covering fraction (C f ) in Table 5. This is calculated from the average of C f (v) in the range of [V out -σ Vout ,V out + σ Vout ] for both the Si II and Si IV multiplet/doublet lines. This choice is because most of the outflowing material (i.e., N ion ) has velocities around the center of the troughs. Therefore, C f represents the covering fraction for the bulk of outflowing material we observe from Si II and Si IV.
Finally, we calculate N ion for Si II and Si IV as below (e.g., Savage & Sembach 1991;Edmonds et al. 2011).
whereN ion (v) is the average column density per velocity over the aperture at v, and N ion is the integrated column density over all velocity bins. The derived N(Si II) and N(Si IV) values for each galaxy is presented in Table 5.
For other singlet transitions such as Si III λ 1206, even though we have the absorption troughs observed in most of our galaxies, the corresponding N(Si III) is difficult to measure due mainly to 3 reasons: 1) Si III λ 1206 is close to Lyα, where the damped Lyα absorption from Milky Way (and sometimes also from the observed galaxy) can cause significant contamination and/or make the continuum hard to determine; 2) Si III λ 1206 is a strong line and usually saturated (see Section 5.3). In this case, the measured N(Si III) from the AOD method will be much smaller than the actual N(Si III) considering the PC effects of outflows; 3) Si III λ 1206 is only a single line. It is not possible to solve the PC equations [i.e., Equations (2) or (3)] and determine robust N(Si III). Overall, we do not present direct measurements of N(Si III) from the spectra. We instead constrain N(Si III) from the photoionization models discussed in Section 5.3.
In Figure 6, we show in the bottom-right panel the resulting N(Si II)(v) for one of our galaxies (J0055-0021). The column density peaks near the line center and declines towards both higher and lower velocities. While the low velocity gas in the outflow will be affected by the accuracy to which the static ISM component is removed, the decline in column densities towards higher outflow velocity is robust. The top-right and lower-left panels show that the decline in N(Si II)(v) is driven by a decline in both the optical depth and the covering fraction.
We also see that the Si II absorption troughs with smaller f (i.e., Si II λ1304) have narrow widths (see top-left panel), which is consistent with Wang et al. (2020). This is as expected from the curve of growth, i.e., while N(Si II)(v) drops in the line wings, stronger lines are more optically thick and, therefore, show wider troughs (see the third panel of Figure  6). We further compare the V out and FWHM out values from Si II λ1260 and 1304 for all our galaxies in the bottom two panels of Figure 4. We find troughs from Si II λ1260 indeed commonly have larger FWHM out than that of Si II λ1304, while their V out are quite consistent. These results imply that the FWHM out parameter is quite sensitive to the distribution of column densities of the outflows (unlike V out ), and there is a relatively narrow range around the characteristic outflow velocity where the bulk of the outflowing gas resides. This gas is traced most clearly by the least optically-thick lines.

Dust Depletion of Gas-Phase Metals
To estimate the total hydrogen column density (N H ) of outflows from the measured N(Si II) and N(Si IV), we need first to determine the amount of depletion of gas-phase silicon onto dust. As discussed in Savage & Sembach (1996); Jenk- . Zoom-in on the spectral regions around S II λλλ1250, 1253, 1259 from the stacked spectra of 31 galaxies that cover all these lines. The normalized spectra are shown in black while the errors are shown in gray at the bottom. The trough from 1259 is commonly blended with Si II λ1260. Therefore, we exclude it when calculating the column density of S II (see Section 5.2).
ins (2009), sulfur has been treated as a standard for an element with very little depletion. Therefore, in this subsection, we compare the measured silicon column densities to those of sulfur to estimate the dust depletion of silicon.
For > half of our galaxies, the HST/COS spectra cover weak lines from S II λλλ1250, 1253, 1259. To measure these weak features in our COS sample, we create a high S/N stacked spectrum that covers these lines (e.g., Alexandroff et al. 2015). A total of 31 galaxies are coadded following the methodology in Thomas (2019). After de-redshifting and re-sampling wavelengths for all galaxies into the same wavelength array, we calculate SNR 2 weighted normalized flux for each bin. We also conduct a 3-sigma clipping in the weighting process, i.e., for one bin, we iteratively remove galaxies with flux that is out of the range [med -3σ, med + 3σ], where med and σ are the median and standard deviation of all flux at this bin in the coadding sample, respectively. This process helps remove outliers in each wavelength bin. The S II region from the stacked spectra are shown in Figure 7.
We follow the same methodology in Section 5.1 to measure the column densities of S II [i.e., N(S II) obs ]. For S II λ1259, since it is too close to Si II λ1260, we exclude it from the calculation. In Table 2, we compare the measured N(S II) from the stacked spectra (4th columns) to the predicted values (5th columns) by scaling N(Si II) as follows: where S/Si is the abundance ratio, and we assume solar abundance; N(Si II) is the mean value of measured N(Si II) for all stacked galaxies (Section 5.1); and ICF(S II) and ICF(Si II) Note. -(1) The number of galaxies stacked.
(2) The predicted column densities are scaled from N(Si II) or N(Si IV) assuming solar ratios of S/Si and typical ionization corrections (see Equation (6)).
are the ionization corrections for S II and Si II, respectively. For typical parameters of our galaxies (see Section 5.3), we get ICF(S II)/ICF(Si II) ∼ 0.6 from CLOUDY models (Ferland et al. 2017). This is consistent with Hernandez et al. (2020) for similar star-forming galaxies. Our choice of a solar abundance ratio for sulfur-to-silicon is a good one because both are alpha elements (created in core-collapse supernovae, Steidel et al. 2016;Kobayashi et al. 2020).
From Table 2, we show that N(S II) obs and N(S II) pred from our stacked spectra are consistent within 1σ. Similarly, we conduct another stack which includes all galaxies covering S IV λ1062 region, and find N(S IV) obs is also consistent with N(S IV) pred within 1σ (see the second part of Table 2). In this case, N(S IV) pred is scaled from N(Si IV) given a solar S/Si ratio and ICF(S IV)/ICF(Si IV) ∼ 0.8 . This suggests that the dust depletion of silicon is similar to sulfur in our galaxies, implying both are primarily in the gas-phase. Qualitativelysimilar depletion patterns have been observed in the Milky Way halo (Savage & Sembach 1996;Jenkins 2009) and in shocked regions in the disk (Welty et al. 2002).
Overall, we conclude that silicon is mostly undepleted into the dust for galaxies in our sample. Therefore, it is viable to estimate N H of outflows from the measured N(Si II) and N(Si IV), which is discussed next.

Methodology
Given the low dust depletion of silicon and the measured N(Si II) and N(Si IV), we can estimate two important properties of the observed outflows, i.e., the total silicon column density (N Si ) and total hydrogen column density (N H ). We do so by running a variety of grid models adopting CLOUDY [version c17.01, (Ferland et al. 2017)]. The fixed parameters are: 1) Starburst99 (Leitherer et al. 1999) models as the input spectral energy distribution (SEDs), where we assume Geneva tracks with high mass loss and constant star-forming history (SFH) with age = 5 Myr. We also assume a standard Kroupa IMF (with slopes of 1.3 and 2.3 in mass ranges of 0.1 -0.5 and 0.5 -100 M , respectively, see Kroupa 2001); 2) an electron number density (n e ) = 10 cm −3 for typical star-burst galaxies; and 3) for each galaxy, we adopt the GASS10 abundance (Grevesse et al. 2010) scaled by the measured O/H values discussed in Section 4.4. The conversion from N Si to N H assumes that the outflow has the same metallicity as measured from the nebular emission lines. 3 We have tested models for n e = 100 cm −3 , and find that the resulting N Si and N H only have minor changes (< 3 -5%). In our grid models, the two varied parameters are: 1) the logarithm of ionization parameter, log(U H ), in the range between -4.0 and 0.0 with a step size of 0.05 dex, and 2) the logarithm of N H in the range between 18.0 and 23.0 [log(cm −2 )] with a step size of 0.02 dex.
For each velocity bin, we solve the best-fit model, i.e., a set of U H (v) and N H (v) values, by comparing the CLOUDY model predicted N(Si II)(v) and N(Si IV)(v) to the measured ones at this velocity (see Section 5.1). This is done through χ 2 -minimizations of the difference between the model predicted and the measured column densities (e.g., Borguet et al. 2012;Xu et al. 2019). Then for the whole outflow, we integrate N H (v) over all velocity bins to get N H . Hereafter, we use N ion (v) to represent the column density per velocity at v [in unit of cm −2 /(km s −1 )], while N ion stands for the integrated column density over all velocity bins in the outflow trough [in unit of cm −2 ]. The best-fit N H values are listed in Table 5, and its distribution is shown in Figure 8, where we find that the mean N H for the observed outflows is 10 20.70 cm −2 . Based on this best-fit model, we can then predict the column densities for other ions , e.g., N(Si III)(v) and N(H I)(v). We also show their integrated values in Table 5. We find that, in all galaxies, Si III is the dominant ion of silicon for the observed outflows, so we estimate the total column density of silicon as N Si = N(Si II) + N(Si III) + N(Si IV).
Furthermore, the mean value of N(Si III) is 10 15.53 cm −2 . For an average FWHMout ∼ 300 km s −1 (see Figure 3), and assuming the trough is flat, we can calculate the average optical depth of Si III by solving Equation (5) as: where we have f Si III = 1.67 and λ Si III = 1206.51 Å. This leads to τ (Si III) ∼ 60. Combined the fact that the observed troughs of Si III are usually non-black (i.e., I > 0), this suggests that Si III troughs of our galaxies are strongly non-black saturated and affected by PC effects. This is consistent with our claims in Section 5.1.
We also see that only a few outflows are consistent with a unit covering factor (with a median value 0.7). The out-flows are therefore somewhat patchy. We see no evidence for any systematic difference between the values of the covering fractions for Si II and Si IV, implying that the two ions likely trace similar regions in the outflow (see also Chisholm et al. 2016a).

The Neutral Phase of the Outflow
For our galaxies, the mean N(H I) from CLOUDY models is 10 18.28 cm −2 , which is only ∼ 0.4% of the mean N H value (10 20.70 cm −2 ). This suggests that our outflows detected in UV absorption lines are mostly ionized gas, and the neutral gas is only a minor part of the outflows in our galaxies.
This also has interesting implications for the origin of the gas traced by low ionization transitions of Si II and C II. These are sometimes used as proxies for neutral hydrogen (e.g., Jones et al. 2013;Alexandroff et al. 2015;Gazagnes et al. 2018). However, the ratios of the observed column densities in Si II vs. those derived via CLOUDY in H I are much too large for the Si II to arise primarily in the H I gas.
More quantitatively, we find that typically only 1 to 10% of the observed Si II comes from the H I phase, so Si II is a poor proxy for the neutral gas. This means that the O I λ1302 line should be used instead, since the nearly identical ionization potentials of O I and H I ensure that the two species arise in the same phase.

Scaling Relationships
In this subsection, we present various empirical scaling relationships between the outflow and galaxy properties. These correlations for low-redshift galaxies (z 0.4) have already been discussed in previous publications in the literature (e.g., Heckman et al. 2000;Martin 2005;Rupke et al. 2005;Heckman et al. 2015;Chisholm et al. 2016a). However, we have conducted more robust analyses for the CLASSY sample. Namely: 1) For each galaxy, our outflow velocity and FWHM are determined from up to 10 lines where both lowand high-ionization transitions are considered (see Section 4.2). In contrast, most previous publications usually have access to only 1 -2 low-ionization lines (e.g., Na I D studied in Heckman et al. 2000;Martin 2005;Rupke et al. 2005). 2) We solve for the outflow's column density given a partial covering (PC) model (see Sections 5.1 and 5.3). In contrast, previous publications are usually aware of the PC property of outflows, but did not include it in their analyses (except in a few cases, e.g., Chisholm et al. 2016b;. This was usually due to the inability to solve the PC equations without doublet or multiplet transitions being well-detected, 3) We have used CLOUDY models to calculate the total H column densities based on the measured ionic column densities. These ionization corrections were not possible in the prior work based on only low-ionization lines, e.g., Na I D and Mg II.
Finally, the CLASSY sample was specifically designed to span maximum ranges in the fundamental galaxy properties of stellar mass, star-formation rate, and metallicity. This makes it ideal for exploring how the outflow properties correlate with these other properties. Overall, it is clearly important to revisit these scaling relationships.
As discussed in Section 3, for all figures in this subsection, we have added 5 LBA galaxies from Heckman et al. (2015) that have HST/COS observations but were not part of the CLASSY sample. These LBA galaxies satisfy all selection criteria of the CLASSY ones, but they can provide additional starbursts with relatively larger values of SFR and M . To be consistent, we follow the same methodology as discussed in Sections 4 and 5 to measure the required quantities for these 5 galaxies. In all figures, galaxies from CLASSY and Heckman et al. (2015) are in red and blue colors, respectively.

Outflow Velocity and FWHM
We begin by examining the correlations of V out and FWHM out (Section 4.2) with the principal properties of the galaxy, including the circular velocity (V cir ), stellar mass (M , in units of M ) and star formation rate (SFR, in units of M /yr). We show these correlations in logarithm scale in Figure 9. The corresponding error bars are shown as crosses. Galaxies with non-detections of outflows are shown as hollow circles, and we arbitrarily set their values to V out = 40 km s −1 and FWHM out = 100 km s −1 (in order to show them in the figures). Note that the non-detections are preferentially in galaxies with low stellar masses and SFRs (for such galaxies, the fraction of detected outflows is smaller).
For each panel, the results of Pearson correlation coefficients (PCC) are shown at the bottom-right corner. The linear fit of the y to x values is shown as the orange dashed line, and the fitted slope and intercept are shown in the top-left corner. We have also conducted the Kendall τ test to assess the statistical significance of each correlation. Both the Kendall and PCC coefficients are listed in Table 4. Objects with UV radii r 50 > 1.5 are labeled as hollow triangles. As discussed in Section 4.4, most of the UV light in these galaxies lies outside of the COS aperture, so their observed absorption outflows may not be representative of the global properties of the galaxy. This is further confirmed by the fact that these galaxies do not lie upon the locus defined by other more centrally concentrated galaxies in Figure 9. Therefore, we exclude galaxies with UV radii r 50 > 1.5 when calculating the correlation coefficients.
As shown in the left column of Figure 9, there are statistically significant correlations between V out (or FWHM out ) with M (and hence with V cir ). It is noteworthy that the slopes of the relationships between V out and FWHM out with V cir are sub-linear, meaning that their ratios decrease with increasing V cir . We will discuss this further in Section 6 below. In the Figure 8. Two distributions for galaxies from our combined sample of CLASSY and Heckman et al. (2015). Left: Distribution for the mean covering fraction of outflows. For each object, this is calculated as the average of the covering fractions for Si II and Si IV around the measured outflow velocity, i.e., Vout ± σVout (see Section 5.1). Right: Distributions for the total column density of outflows (NH, in units of cm −2 ), which are derived from the photoionization models of CLOUDY (see Section 5.3). Note that for both panels, galaxies labelled as "No outflow" are not included. right column of Figure 9, we see strong positive correlations between V out (or FWHM out ) with SFR. The strengths of the correlations with SFR and M are very similar. This is in part because SFR and M are themselves well-correlated in this sample (Paper I). In contrast, the sample studied in Heckman et al. (2015) included galaxies observed by Far Ultraviolet Spectroscopic Explorer (FUSE) that had significantly lower values of SFR/M (typically 10 −10 to 10 −9 M ) compared to CLASSY. With the inclusion of these galaxies, they found a stronger correlation of V out with SFR than with M .
We have also noticed that the correlations with FWHM out are slightly stronger than the ones with V out in Figure 9. As we explained in Section 5.1, unlike V out , FWHM out is sensitive to both the bulk kinematics of the outflow and to the distribution of column density (as reflected in the extent of the broad and shallow wings of the outflow profiles). These wings represent lower column densities in outflows and are more sensitive to the variations of column densities (see Figure 6). Therefore, differences in galaxy properties could affect FWHM out somewhat differently than V out .
In Figure 10, we also test the correlations between the normalized outflow velocity (V out /V cir ) and normalized measures of the star formation rate (SFR/area and SFR/M ). Previous studies found positive correlations for them (e.g., Martin et al. 2012;Heckman et al. 2015). Our galaxies in these figures show large scatter, and the correlations are weak or insignificant. As noted above, the main difference between our sample and the Heckman et al. (2015) sample is that we do not have galaxies with the relatively low values of SFR/area and SFR/M represented by their galaxies with FUSE data. In a future paper, we will re-analyze the FUSE data in the same way as we have done for the current sample and then revisit these potential correlations.
Similar results relating V out to SFR, M , and V cir have been found previously in low-z starbursts (Martin 2005;Rupke et al. 2005;Chisholm et al. 2015;2016a;Heckman & Borthakur 2016). A direct comparison to our results is not straightforward, largely because few of the studies defined V out the same way we have done. In some cases, a "maximum" velocity was used rather than a line centroid (e.g., Rupke et al. 2005;Heckman & Borthakur 2016). In other cases, the entire absorption feature was treated as a single component (e.g., Heckman et al. 2015;Chisholm et al. 2015;2016a). The only study that used something similar to our double-Gaussian approach was Martin (2005). Another difference is that Martin et al. (2012) and Rupke et al. (2005) used the Na I D doublet to probe the outflows. This traces a dusty H I component in the outflow, while the UV lines used in CLASSY, Chisholm et al. (2015), Chisholm et al. (2016a), and Heckman & Borthakur (2016) trace warm ionized gas.
With these caveats in mind, we summarize the results from these papers and compare them to CLASSY in Table 3. Given the differing definitions of V out , we only list the loglog slopes (and not their normalizations). All the studies are rather consistent. The main advantage of the CLASSY sample is an improved sample size at low values of SFR, V cir , and M .

Outflow Rates: General Considerations
Before using these data to estimate outflow rates in the galaxies, it is useful to briefly consider the general methodology and resulting uncertainties. We will use the mass outflow Figure 9. The log of the outflow velocity (Vout) and Full-Width-Half-Maximum (FWHMout) plotted as a function of the basic properties of the star-forming galaxies. From top-left to bottom-right are log(Vout) vs. circular velocity, star formation rate, and log(FWHMout) vs. circular velocity, star formation rate, respectively. Galaxies with HST/COS observations from CLASSY and Heckman et al. (2015) are shown in red and blue, respectively. The corresponding error bars are shown as crosses. The gray solid lines in the left panels indicate Vout= Vcir and Vout= 10Vcir (and the same for FWHMout). Galaxies with non-detections of outflows are shown as hollow circles. We set their Vout = 40 km s −1 and log(FWHMout) = 100 km s −1 (only to include them in the figures). Galaxies with large UV sizes compared to COS aperture are labeled as hollow gray triangles. For these galaxies, the COS aperture covers only part of the starburst and the data may reflect local rather than global outflow properties. These galaxies' spectra will also be more affected by vignetting. We therefore exclude them when calculating the correlation coefficients. The results for Pearson correlation coefficients (PCC) are shown at the bottom-right corner of each panel. The linear fit of the y to x values is shown as the orange dashed line, and the fitted slope and intercept are shown in the top-left corner. See discussion in Section 5.4.1. For the top-right panel, we also show the best-fit wind-blown bubble model in green assuming momentum-conserving, which matches the data well [see Section 6.2 and Equation (18)]. rate as a specific example, but these general considerations will apply to all the outflow rates discussed later.
Simple dimensional analysis tells us that the average mass outflow rate (Ṁ out ) will just be the mass of outflowing gas (M out ) within a radius R out divided by the time it takes the flow to traverse this distance, i.e., R out /V out . This leads tȯ M out = M out V out /R out , where V out is the outflow velocity. The strength of the observed absorption-lines depends on the column density of the outflowing gas (particles per unit area).
Let us consider a few idealized cases. In the simplest case of an expanding thin shell with mean velocity V shell , radius R shell , and solid angle Ω we have: where µ is the mean mass per proton (∼ 1.4) and m p is the proton mass. The difficulty in using this equation is that we do not know the value of R shell . We will return to consideration of shells in the context of a wind-blown bubble model in Section 6 below. Figure 10. Comparisons of normalized outflow velocity (Vout/Vcir) with normalized star formation rate (SFR/M and SFR/area). The captions and labels are the same as Figure 9. Due to the small dynamic range of Vout/Vcir in our sample, the correlations in these two panels are weak (L) to insignificant (R). See discussion in Section 5.4.1. Figure 11. Correlations related to the silicon mass flow rate (ṀSi,out) and total silicon column density (NSi). Labels and captions are the same as Figure 9. Left: We compareṀSi,out from the outflow to the amount that is provided from the starburst. The solid lines from left to right represent where y = 10, 1, 0.1 x, separately. Right: We compare NSi to the dust extinction calculated from fitting the stellar continuum (Section 4.4). We also show the prediction (in black line) from Draine (2011) assuming a standard Milky Way dust/metals ratio and the observed median dust extinction in outflow ∼ 16% of the total dust measured. The right y-axis represents the dust optical depth at 1200 Å assuming Reddy et al. For now, we consider cases of continuous mass-conserving outflows with different radial profiles of density and velocity. The simplest such case is an outflow with constant velocity V out and a density profile n(R) = n 0 (R/R 0 ) −2 , where R 0 is the radius at which the outflow begins. This simple case is broadly consistent with both numerical simulations (e.g., Schneider et al. 2020) and analytic models (e.g., Fielding & Bryan 2022) of multi-phase galactic outflows. In this situation (hereafter, case 1), it is straightforward to show thaṫ However, based on the analysis of the UV emission-line properties of outflows (e.g., Wang et al. 2020;Burchett et al. 2021), we also consider a shallower radial density profile n(R) = n 0 (R/R 0 ) −1 , which implies V (R) = V max (R/R 0 ) −1 for mass-conservation (and also means that the outflow's maximum velocity is at R 0 ). In this case (hereafter, case 2): where R max is the maximum radius which the outflow reaches. -We compare the slopes of the scaling relationships with other known samples of star-forming/starburst galaxies in the literature. The log-log slopes between Vout and SFR, Vcir, and M are shown in the 5th, 6th, and 7th columns, respectively.
(2). Definitions for Vout. In this paper, we derive Vout from the double-Gaussian fitted results (see Section 4), which is similar to Martin (2005). Other publications adopted either single-Gaussian fitting or took the maximum velocity (or v90) of the troughs as Vout.
(3). From HST/COS spectra, the major rest-frame Far-UV lines adopted to estimate Vout are: O I, C II, Si II multiplet, Si III, Si IV (used in this paper, see Section 4, and Chisholm et al. 2016a;Heckman & Borthakur 2016). For works that also considered FUSE data (Heckman & Borthakur 2016), additional Far-UV lines are adopted, including C III λ977, C II λ1036, and N II λ1084.
(4). These did not exclude the cases in which the COS aperture did not cover at least 50% of the starburst FUV continuum (see Section 4.4).
The key point here is that for reasonable choices for the relevant parameters, the values ofṀ out will be the same to better than a factor of two in two radial density laws described above. More explicitly, if we define V max ∼ V out + 0.5 FWHM out , then Figure 5 implies V max ∼ 2 V out . The mass outflow rates will be exactly the same in case 1 and case 2 above if ln(R max /R 0 ) = 2 (i.e. R max = 7.3 R 0 ). Since this involves the log of the ratio, its value depends only weakly on R max /R 0 .

Metal Mass Outflow Rates and Dust Extinction
With these considerations in mind, we calculate the metal mass outflow rates of silicon (Ṁ Si,out ) for galaxies in our sample. One advantage ofṀ Si,out over the total (i.e., hydrogen) mass outflow rate (Ṁ out ) is that we do not need to know the metallicity of the outflows, for which we don't have direct measurements.
Let us consider the simple cases of an outflow with a constant velocity and a R −2 density profile [see Equation (9) where dN Si /dv = N Si (v) is the silicon column density per velocity, and the integration is for the velocity range of the observed outflow trough. We have three unknowns (i.e., Ω, N Si , and R 0 ) required to calculateṀ Si,out . For Ω, as shown in Section 4.2, in the CLASSY sample, ∼ 85% of galaxies are identified as "hosting outflows". Given this high detection rate and the possible existence of outflows with directions in the sky-plane (which are undetectable), we take Ω = 4π. The value for N Si is derived from the CLOUDY models discussed in Section 5.3. For R 0 , previous studies either assume a fiducial radius (e.g., 1 -5 kpc in Martin 2005;Rupke et al. 2005;Martin et al. 2012), or relate it to the starburst radius (e.g., Chevalier & Clegg 1985;Heckman et al. 2015). We take the second method and assume R 0 = 2 × r 50 (the outflow begins at twice the half-light radius of the starburst). This choice is consistent with Heckman et al. (2015), and corresponds to an outflow that begins at a radius enclosing 90% of the starburst (for an exponential disk model). It is also consistent with the radius at which the outflow of ionized gas is observed to begin in M 82 (Shopbell & Bland-Hawthorn 1998). As noted above, the outflow rates scale linearly with the adopted radius, so this choice is important.
In the left panel of Figure 11, we compareṀ Si,out with the rate at which silicon is created and injected by SNe in the starburst, i.e.,Ṁ Si,* (both in units of M /yr). We approximatė M Si,* = 1×10 −3 SFR assuming Starburst99 models for the Si yield (Leitherer et al. 1999;Heckman et al. 2015). While there is a good correlation, the median ratio ofṀ Si,out anḋ M Si,* is only ∼ 0.2. This implies inefficient incorporation of the Si ejected by SNe into the warm ionized outflow that we observe. We will discuss this further in Section 6.4 below.
In the right panel, we show the relationship between N Si and E(B − V ). The latter is the dust extinction derived from stellar continuum fits to the FUV spectra for each galaxy (see Section 4.4). There is a positive correlation, but with considerable scatter. A correlation is not surprising since both N Si and dust extinction are related to the total column density of metals. The difference is that N Si only measures metals in the outflow, while the extinction includes dust in both the outflow and the static ISM. Furthermore, the scatter is also not surprising. For one thing, N Si is derived from partial covering fraction models, while E(B-V) is for a screen that covers the starburst.
To evaluate this further, we can estimate the maximum amount of dust extinction in the outflow using N Si , assuming a standard Milky Way (MW) dust/metals ratio (Draine 2011) and assuming the solar ratio of silicon to metals mass. Figure 12. Correlations related to total mass flow rate (Ṁout) and mass loading factor (Ṁout/SFR). Labels and captions are the same as Figure  9. The top left panel shows that the mass outflow rate has sub-linear dependence on the SFR. The top right and lower left panels show how the ratio ofṀout/SFR scales with the galaxy circular velocity and with the outflow velocity. In both cases, the ratio increases rapidly at lower velocities (slopes ∼ -1.6 to -1.7). See further discussion in Section 5.4.4. This then implies that the total dust extinction is related to N Si as (see Draine 2011): In our sample, we then check the distribution of N Si /E(B-V) tot /(1.9 × 10 17 ), and the resulting median is ∼ 0.16. This suggests that the dust in our observed outflows only represents ∼ 16% of E(B-V) tot . Therefore, most of the dust responsible for the observed extinction must reside in the static ISM. This relationship [Equation (12) with multiplication factor of 0.16 on the right side] is shown as the black line in Figure 11. On the right y-axis, we also show the implied upper limits on the dust optical depths in the outflow (τ 1200 ) at 1200 Å assuming the Reddy et al. (2015) extinction law and Equation (12). This leads to a typical value of < τ 1200 > ∼ 0.3 due to dust in the outflow.

Total Mass Outflow Rates and Mass Loading Factors
Similarly to the calculations ofṀ Si,out [see Equation (11)], we can derive the total mass outflow rate (Ṁ out ) as: where v is the outflow velcity, dN H /dv = N H (v) is the total hydrogen column density per velocity and has been derived in Section 5.3. We use the same values described above for R 0 and Ω. In Figure 12, we present the correlations related tȯ M out .
In the first panel, we show a strong correlation betweeṅ M out and SFR, which has been found in previous publications for low-redshift galaxies (e.g., Martin 2005;Rupke et al. 2005;Heckman et al. 2015). For our combined sample,Ṁ out ranges from 0.3 -40 M /yr, and the best linear fit yields thaṫ M out = 3.8 × SFR 0.41 (a sub-linear slope).
In the next two panels, we show the correlations between the mass loading factor (i.e.,Ṁ out /SFR) and V cir (and V out ). Both figures show strong inverse correlations, with the correlation betweenṀ out /SFR and V cir being stronger. The best fit slopes for both figures are ∼ −1.6 to −1.7, intermediate between the expectations for a so-called "momentum-driven" and "energy-driven" outflows. We will discuss this in Section 6. It is striking that the mass-loading factors reach values of ∼ 10 in the lowest mass galaxies.

Momentum and Kinetic Energy Outflow Rates
Given the calculatedṀ out , we can derive the momentum and kinetic energy outflow rates as: (14) The integration over v inṗ out andĖ out is necessary since both functions are strongly dependent on v. We can then compareṗ out andĖ out with the momentum and kinetic energy supplied from starburst, i.e.,ṗ andĖ . As discussed in Heckman et al. (2015),ṗ is a combination of a hot wind fluid driven by thermalized ejecta of massive stars (Chevalier & Clegg 1985) and radiation pressure (Murray et al. 2005). In Starburst99 models (Leitherer et al. 1999), this leads tȯ p = 4.6×10 33 SFR dynes, where SFR is in units of M /yr. Similarly,Ė = 4.3×10 41 SFR erg/s. For each galaxy, these values are listed in Table 5.
We present the comparisons in the first two panels of Figure 13, where we find strong positive correlations for both. We draw several solid black lines to represent where Y = 10, 1, 0.1 X in the left panel and Y= 1, 0.1, and 0.01 X in the right panel. For most of galaxies in our combined sample, we haveṗ out = 10 -100%ṗ , andĖ out = 1 -20%Ė . These ranges are similar to the ones reported in ?, but both several times smaller than the typical values derived by Heckman et al. (2015) based on a much simpler analysis.
In both cases, the slopes are sub-linear, implying that as the SFR increases, less of the available momentum and kinetic energy supplied by stars is being carried by the warm ionized phase of the outflow that we are probing. This trend is stronger for the momentum than for the kinetic energy (see more discussion in Section 6).
We see that there are no cases in whichĖ out >Ė . However, there are three galaxies in whichṗ out >ṗ . This can be understood in the situation in which the wind fluid that accelerates the gas we see has been substantially mass-loaded but conserves kinetic energy. We will return to this in Section 6.

DISCUSSION
We begin the discussion by summarizing the most widely used simple model for galactic winds. We then describe several straightforward analytical models for the outflows. This is then followed by a discussion of the scaling relations often adopted for galactic winds in semi-analytic models and numerical simulations of cosmic volumes and we contrast them to our results. Finally, we compare our results to a recent analytic models and high-resolution numerical simulations of multi-phase galactic winds.

Theoretical Background
A very influential model for galactic winds driven by a population of massive stars was that of Chevalier & Clegg (1985). It was simplified, with an assumption of spherical symmetry and neglected gravity. It was also a single-phase model treating only the thermalized ejecta of massive stars (supernovae and stellar winds). As such, it served as a important foundation for more elaborate models to come. The model predicts a very hot region of gas within the starburst. This gas then passes through a sonic point at the starburst radius and then transitions into a supersonic galactic wind. The inputs to the model are the mass and kinetic energy injection rates from the massive stars (Ṁ andĖ ) and the radius of the spherical starburst. For the purposes of this paper, the most important outputs of the model are the terminal velocity of this wind, and the outflow rates (of mass, momentum, and kinetic energy, all of which are conserved with radius). These quantities are given as: In the above equations β accounts for ambient gas that is mixed into the stellar ejecta (so that the total outflow rate in the wind fluid is a factor β larger than the rate at which stellar ejecta are created). 4 The term α accounts for the effects of radiative losses that drain away the energy carried by the wind fluid. Thus, β ≥ 1 and α ≤ 1.
It is crucial to emphasize that the outflowing gas described above is not the gas we observe in absorption. In the rest of the discussion we will make this distinction by referring to the former as the wind fluid and the latter as the warm outflow. In this simple standard model, the warm outflow traces preexisting gas clouds that are being accelerated via momentum and/or kinetic energy transferred from the wind fluid to the clouds 5 . The initial idea (Chevalier & Clegg 1985) was that the clouds were accelerated by the ram pressure of the wind fluid. A challenge has been understanding how clouds survive being shocked by the wind fluid long enough to be accelerated to the observed velocities in the warm outflows (Heckman & Thompson 2017).
Note that radiation pressure from the starburst coupled to dust in the clouds can also transfer momentum to the clouds (Murray et al. 2005). However, for the choices α = β = 1, the 5 There is an alternative class of models in which the wind fluid suffers strong radiative cooling and the warm outflow forms directly out of the wind fluid (Thompson et al. 2016;Schneider et al. 2018;Lochhaas et al. 2021). We will compare these models to the data in Section 6.6 below.
wind fluid's ram pressure is about four times larger than the radiation pressure (L bol /c, where L bol is the bolometric luminosity). We have seen in Section 5.4.3 above that any dust in the outflow is optically-thin in the far-UV (τ 1200 ∼ 0.3). This means that only ∼ 25% of the momentum provided by UV radiation can be transferred to dust in the outflow. Thus, P ram /P rad > 16, and radiation pressure is most likely negligible in driving these outflows. One relatively new and promising idea is that the momentum transfer from the wind fluid to the clouds occurs in turbulent mixing layers at the interface between the cool cloud and the wind fluid, and (more generally) that mass and momentum can be exchanged in both directions between the clouds and the wind fluid (e.g., Gronke & Oh 2020;Fielding & Bryan 2022). We will explore such models in Section 6.4 below.  We can now examine some of the scaling relations we have measured in the context of the simple standard model described above of an outflow driven by a fast-moving wind fluid. Observations of the hard X-ray emission in the M 82 starburst favor the choice α ∼ β ∼ 1 (Strickland & Heckman 2009). In this case, V wind ∼ 2700 km s −1 . This is over an order-of-magnitude higher than the warm outflow velocities we see, but again, we are not directly observing the wind fluid itself.
One key prediction of this model is that the kinetic energy carried by the warm outflow cannot exceed the kinetic energy supplied by the starburst. Our results in Figure 12 are fully consistent with this, as the typical ratios ofĖ out /Ė are 1 to 20 %. We note that this does not imply a small value for α (which pertains to the wind fluid rather than the warm outflow). Indeed the models and simulations discussed in Sections 6.3 and 6.4 below find that the great majority of the kinetic energy is carried by the wind fluid rather than the warm outflow.
The momentum flux in the warm outflow reflects the momentum transferred to it from the wind fluid. We see in Figure 13 that the measured momentum fluxes are generally similar to (but smaller than) the predicted momentum flux in the wind fluid, with a median ratio of ∼ 30%. This suggests an efficient transfer of momentum. Interestingly, there are several cases in whichṗ out >ṗ . This is not in conflict with the model of an outflow driven by the wind fluid. That is, in the case where α × β > 1, the momentum flux in the wind fluid can exceed the momentum flux from the starburst by a factor (α × β) 1/2 (see Equation set (16)). This would correspond to a case in which radiative losses in the wind fluid are negligible (α ∼ 1) but the wind fluid is significantly contaminated by ambient gas mixed into it (β >> 1). This would of course require an efficient transfer of this momentum from the wind fluid to the clouds in the warm outflow.

Comparison to a Simple Analytic Models of Momentum-Driven Outflows
We first compare our observations with a model of a population of clouds acted upon by a combination of outward momentum from the starburst and the inward force of gravity, as discussed in Heckman et al. (2015). They have derived a critical momentum flux required to drive the outflow (see their Equation 3): The resultingṗ crit values are listed in Table 5. In the last panel of Figure 13, we compare the normalized outflow velocity, i.e., V out /V cir , versus the normalized outflow momentum flux, i.e.,ṗ out /ṗ crit . The two black vertical lines split the figure into three different regimes for outflows: (1)ṗ out /ṗ crit < 1, where no outflow is expected to be driven due to the lack of momentum inputs. This is consistent with the fact that no outflows from our combined sample fall in this region. (2) 1 <ṗ out /ṗ crit < 10, where we expect relatively weak outflows are driven in this regime. A few of our observed outflows fall into this region.
(3)ṗ out /ṗ crit > 10, where relatively strong outflows should be driven. We find most of our observed outflows fall into this last region ofṗ out /ṗ crit . This is as expected since weaker outflows are harder to detect or more easily marked as "no outflows" due to various issues mentioned in Section 4.2.
In in last panel of Figure 13, we also show the expectations from Equation (5) of Heckman et al. (2015), where they predict the maximum velocity of an outflowing cloud (V max ). This velocity corresponds to the radius at which the outward force due to ram and radiation pressure equals the inward force due to gravity (and beyond which the cloud begins to decelerate). In dotted, solid, and dashed curves, we overplot this equation on our data where we assume the observed outflow velocity V out = 100%, 50%, and 20% of V max , respectively. We find that most of our confirmed outflows are consistent with the V out = 50% V max curve. This can be explained as a "natural choice", i.e., V max = V out + 0.5 FWHM out ∼ 2V out (see Section 5.4.2). Unlike Heckman et al. (2015), we do not probe the regime of lowṗ out /ṗ crit very well, and so do not see its correlation with V out /V cir . In a future paper, we will analyze the Heckman et al. (2015) FUSE data using the methods in this paper and revisit this plot.
Next we consider the simple model of wind-blown bubble driven into the ISM/CGM by the momentum flux of the wind fluid. 6 In this model, the gas we see in absorption is sweptup ambient gas at the surface of the expanding bubble. For simplicity, we use the simple spherically-symmetric model described in Dyson (1989) in which the ISM is treated as having a uniform density and gravity is neglected. The latter is supported by the large values ofṗ out /ṗ crit seen in Figure 13 above.
In this case, the predicted radius and expansion velocity are given by: Hereṗ 35 is the momentum flux in the wind in units of 10 35 dynes, SFR 10 is the star formation rate in units of 10 M yr −1 ) (the median value for the outflow sample), n o is the H density (cm −3 ) in the region into which the bubble expands, and t 7 is the time since the expansion began in units of 10 7 years, similar to the ages derived from fits of SB99 models to the COS data (Senchyna et al. 2022, in preparation).
It is intriguing that the predicted dependence of V s ∝ṗ 1/4 35 (i.e., V s ∝ SFR 1/4 ) is very similar to the measured slope between V out and SFR (∼ 0.24, see Figure 9). In the second panel of Figure 9, we show the best-fit relation with n 0 t 2 7 6 We have also considered the case in which the bubble expansion is driven by the kinetic energy of the wind fluid. This model is a poorer fit the data, so we omit a discussion here. treated as a free parameter (green dashed line). We find a best-fit value of n 0 t 2 7 = 2.8 ×10 −2 cm −3 yr 2 . Assuming t 7 ∼ 1, this implies low-density gas (i.e., 2.8 ×10 −2 cm −3 ), which should be located well outside the starburst region.
We can go further with this model. The predicted column density through the bubble wall is given by n 0 R s /3. For the above value of n 0 and t 7 , Equation (18) yields N H = 1.5 × 10 20ṗ 1/4 35 cm −2 , which is only a factor of two to three below the typical values we measure (see Figure 8). However, we do not see the predicted decline in N H with decreasingṗ * in the data.
Finally, we can compute the mass-outflow rate in this model (the time averaged rate at which ambient gas has been swept-up by the expanding bubble). Using the equations above, and taking n 0 = 2.8 × 10 −2 cm −3 , we get: These rates are about three times larger than we have estimated, and would lead to a steeper slope in the dependence of M on SFR (0.75) than we observe (0.4, see Figure 12).
We conclude that a simple wind-driven bubble model has some success in matching the data. However, it cannot be a complete model: for an expanding bubble, the observed absorption-lines will only come from the part of the bubble located directly along the line-of-sight to the starburst. This will result in a narrow blueshifted absorption-line with FWHM out << V out . This is inconsistent with the data (see Figures 2, 3, and 5).

Theoretical Scaling Relations in a Cosmological Context
There is a substantial literature in which the properties of galactic winds in a cosmological context are modeled using simple physically-motivated prescriptions (e.g., Somerville & Davé 2015;Naab & Ostriker 2017). These include both semi-analytic models and numerical simulations in which the winds are not modeled ab initio (due to insufficient spatial resolution), but rather, are implemented using various subgrid recipes.
Here we briefly summarize these popular prescriptions and compare them to the data. For the most part, these models assume that there is a linear proportionality between V out and V cir (e.g., Guo et al. 2011;Davé et al. 2013). This is not fully consistent with the results shown in Figure 9, which imply V out ∝ V 0.6 cir . At face value, this would suggest that the outflowing gas is more likely to escape from the low-mass galaxies.
There are several different prescriptions for how the massloading factorṀ out /SFR should vary with V out . Since these models assume a linear relationship between V out and V cir , the mass-loading factor will scale the same with both velocities. In one case the prescription is that outflows all carry the same fraction of the momentum supplied by the starburst. In this case,Ṁ out /SFR ∝ V −1 out ∝ V −1 cir (Oppenheimer & Davé 2008;Dutton et al. 2010). An analogous assumption is that the outflows instead carry a fixed fraction of the kinetic energy supplied by the starburst (Baugh et al. 2005;Somerville et al. 2008), leading toṀ out /SFR ∝ V −2 out ∝ V −2 cir . Finally, some models assume that the mass-loading factor is a constant (e.g., Davé et al. 2013;Ford et al. 2016).
As seen in Figure 12, the observed scaling relations of the mass-loading factor with V out and V cir are both intermediate between these two cases (slope ∼ -1.6 to -1.7). Note that Heckman et al. (2015) found a shallower slope. The main difference is that the CLASSY sample reaches much lower values of V cir (providing better dynamic range).

Comparisons to an Semi-Analytic Model of a Multi-phase Outflow
In this section, we compare our results to a new semianalytic model for multi-phase galactic winds (Fielding & Bryan 2022). This model starts with the Chevalier & Clegg (1985) model for the wind fluid, but adds radiative cooling and gravity. More importantly, the model incorporates the bi-directional exchange of mass, metals, momentum, and energy between the wind fluid and a population of much denser clouds. These transfers take place in turbulent mixing layers arising as the wind fluid flows along the cloud surface. It is the interaction between the wind fluid and the clouds that lead to the warm outflow we observe.
The case considered is based on M 82, with V cir = 150 km s −1 , SFR = 20 M year −1 , and a starburst radius of 300 pc. The free parameters that are explored are what we have called β for the wind fluid, which is varied between 1 and 5, the initial mass-loading factor (Ṁ out /SFR) of the warm phase arising as the wind fluid interacts with ambient clouds, which is varied between 0.2 and 0.5, and the mass of an individual cloud, which varies from 10 2 − 10 6 M .
The outflow velocity in the model depends most strongly on the cloud mass, with the most massive clouds being accelerated to the lowest velocity. Our results in Figure 9 imply that V out will be 260 km −1 based on the SFR in the model and 360 km s −1 based on V cir in the model. These velocities are most consistent with the most massive clouds in the models (10 6 M ), which reach V out = 350 (450) km s −1 at radii of 1 (10) kpc.
The final mass-loading factor for clouds in the model depends on both its initial value, and the subsequent transfer of mass between the wind fluid and clouds as they flow out. Our results in Figure 12 imply that an M82-like galaxy would have a mass-loading factor of ∼ 0.6 based on SFR and ∼ 0.3 based on V cir . This is in the range adopted by the models for the initial mass-loading factor. In the models, the massloading factor in the warm outflow only drops significantly with radius for the case in which β = 1 (uncontaminated wind fluid) and the cloud masses are small (< 10 4 M ). In these cases, the clouds are shredded and incorporated into the wind fluid. These particular models are not consistent with the data.
The momentum flux carried by the clouds in the models shows little radial dependence and ranges from ∼ 30 to 100% ofṗ for β from ∼ 1 to 5, respectively. The β = 1 model is therefore a better fit to our data ( Figure 13). Finally, the kinetic energy flux in the clouds depends most strongly on β, being in the range 20 to 40%Ė for β = 5 vs. only 4 to 10% for β = 1. Our median value of ∼5% (see Figure 13) again favors the β = 1 models.
Next, we can consider the question of the mixing of the metals created and dispersed by SNe and carried by the wind fluid with the material observed in the warm outflow. The Fielding & Bryan (2022) models do not explicitly calculate that quantity, but they do compute the separate radial dependences of the metallicity of the wind fluid and of the clouds for the different models. They assume an initial metallicity of twice solar in the wind fluid and solar in the clouds. For strong metal mixing, the wind and cloud metallicities converge at large radii. This is the case for all cloud masses in the β = 5 model, but only for the lower mass clouds in the β = 1 model. In Section 5.4.3, we have found that the observed outflow rate of Si in the warm gas has a median value of only about 20% of the Si creation rate by SNe (consistent with most of the newly-created Si residing in the wind fluid). This result is at least most qualitatively consistent with the β = 1 model with massive clouds.
These models also make different predictions for the total hydrogen column density (N H ) in clouds along the line of sight (D. Fielding, private communication). For the β = 1 models, N H increases from about 2 ×10 20 to 8 ×10 20 cm −2 as the individual cloud masses increase from 10 2 to 10 6 M . For the β = 5 model the corresponding range is ∼ 1 to 2 ×10 21 cm −2 . For our sample we find a median value of N H = 5 × 10 20 cm −2 , which agrees best with the β = 1 model and massive clouds.
In summary, at least for parameters appropriate to an M82level starburst, our data are most consistent with the models with β = 1 (a fast wind that has not been contaminated) and with relatively massive clouds (∼ 10 5 to 10 6 M ). In the future, it will be interesting to compare our data to new models that have been tuned to cover the ranges of the CLASSY sample in V cir , SFR, and starburst size.

Comparisons to a Numerical Simulation of a
Multi-phase Outflow XU ET AL.
First, we compare our results to the highest-resolution numerical simulations of galactic winds currently available (Schneider et al. 2020). Their simulations are designed to roughly correspond to M 82 in terms of SFR and V cir , but the starburst is modeled as a population of massive star-forming clumps distributed within a radius (R * ) of 1 kpc. The results of the simulations are shown for gas in three temperature regimes: 1) a hot phase, i.e., T > 5 × 10 5 K, 2) an intermediate phase with 5 ×10 5 > T > 2 × 10 4 K, and 3) a cool phase with T < 2 × 10 4 K. We focus here on the results shown for a time of 35 Myr (since they significantly drop the input SFR at later times).
The predicted values for V out for the cool phase rise rapidly with radius to a value of ∼ 500 km s −1 at 4 kpc (4 R * ) and then gradually increase to about 750 km s −1 at 10 kpc. Our results in Figure 9 imply that V out will be 260 km −1 based on their assumed SFR and 360 km s −1 based on their assumed V cir , roughly half the values from the simulations.
The outflow rates presented in the figures in Schneider et al. (2020) refer only a bi-conical region with an opening half-angle of 30 • . They find that integrating over 4π ster leads to total outflow rates that are three times larger. We use this scale factor to compare to our results. This yields a mass-loading factor in the cool phase that peaks at about 10% at radii between 2 and 5 kpc and then drops at larger radii. This is significantly smaller than our results in Figure  12, which imply that an M82-like galaxy would have a massloading factor of 0.6 based on SFR and 0.3 based on V cir .
Similarly, the momentum flux in the cool phase in the simulations is about 10 to 20% ofṗ (falling at large radii). This is a bit smaller than the median value of ∼ 30% in the data ( Figure 13). Finally, the kinetic energy flux in the cool phase in the simulations is in the range of ∼ 2 to 4% ofĖ , compared to a median value of 5% in our data ( Figure 11). As in the analytic models by Fielding & Bryan (2022), the lion's share of the kinetic energy is in the hot phase (i.e., the wind fluid).
Overall, the results in the simulations of Schneider et al. (2020) are in reasonable agreement when compared with the galaxies studied here, but overall produce outflows of the cool gas that are too fast and do not carry enough mass, momentum, and kinetic energy (all by factors of ∼ 2). The simulations most resemble the analytic models of Fielding & Bryan (2022) in the regime of low cloud masses.
There are several other publications that present simulations of multiphase galactic winds and predict various scaling relationships. We briefly discuss two of them below, which are both built on Athena MHD code (Stone et al. 2008): 1). Tanner et al. (2017) adopts 3D hydrodynamical simulations assuming similar models as Equations (15) and (16) to predict the outflow velocities for various silicon ions (Si I -Si XIII). They find velocities increase significantly ( a Figure 14. Correlations between Vout and the SFR surface density for galaxies in our combined sample. Labels and captions are the same as Figure 9. The black solid line represents the predictions from multiphase outflow simulations presented in Kim et al. (2020). The slopes match, but the model velocities are about a factor of 3 -5 smaller than the measured ones (see Section 6.5).
factor of 2) from Si II to Si IV, which is not seen in our galaxies ( Figure 4). They also predict that V out versus SFR correlation has an abrupt flattening above SFR = 5 -20 M /yr defined by the hot wind velocity. We do not see any statistically significant evidence of this flattening in SFR > 10 M /yr in Figure 9. 2). Kim et al. (2020) presents a suite of parsec-resolution numerical simulations for multiphase outflows and shows that V out correlates with the surface density of SFR with a slope ∼ 0.2, which is similar to what we get in our combined sample ( Figure 14). However, the normalization of the outflow velocity from the simulations is too low compared to the data by a factor of ∼ 3 -5. These simulations predict that only about 10% of the metals injected by supernovae are incorporated in the warm ionized phase, with the majority being in the hot phase. This is consistent with our results.
All-in-all, there exist both consistencies and differences between our observations and various simulations. To gain more insights into the physics of galactic outflows and their correlations with galaxy properties, incorporation of observations into the recipes of simulations are necessary, which is one of our main goals in future papers.

Comparisons to Models of Cooling Winds
One interesting idea is that the observed outflowing warm gas is produced directly through radiative cooling of the wind fluid. This can happen when there is enough contamination of the wind fluid (β >> 1) that the radiative cooling time of the poisoned wind fluid is shorter than the outflow time of the fluid (e.g., Wang 1995;Martin et al. 2015;Thompson et al. 2016). Figure 15. Comparisons of predicted outflow velocities from models of cooling winds (Thompson et al. 2016) to our measured outflow velocities. Left: Comparisons between the predicted maximum outflow velocity [see Equation (20)] to the measured maximum outflow velocity. Right: Comparisons between the predicted critical outflow velocity [see Equation (21)] to the measured outflow velocity. For both panels, the orange dashed lines represent the best fitting linear correlations, while the black lines show where y = x. The models substantially over-predict the values of the observed outflow velocities. See discussion in Section 6.6. Thompson et al. (2016) consider two cases. One is a case in which β is sufficiently large that radiative cooling becomes important somewhere in the wind (at a radius >> the starburst radius of R * ). This leads to a predicted maximum outflow velocity given by: Here ξ is the metallicity of the wind fluid (relative to solar), SFR 10 is in units of 10 M /yr, and R * ,0.3 is the starburst radius in units of 0.3 kpc.
They also consider a case in which β is large enough for cooling to occur at R * .This yield a critical outflow velocity given by: In Figure 15, we compare these two velocities to our observed values of V max = V out + 0.5FW HM out and V out , respectively. We see that in both cases the predicted velocities are significantly larger than the observed values (by typical factors of 3 to 10). The discrepancies are particularly large for the slower outflows. In simple terms, these slow outflows require such large values of β that the outflows would be stillborn inside the starburst. Also, these models have no natural explanation for the correlation between outflow velocity and galaxy circular velocity (see Figure 9).
The cooling wind models was explored in more detail in high-resolution numerical simulations of a multi-phase starburst-driven wind modeled on the prototype of M 82 (Schneider et al. 2018). The simulations produce outflow ve- Figure 16. Comparisons between the measured momentum flux of outflows (ṗout) to the predicted maximum wind momentum flux (ṗmax) from Lochhaas et al. (2021). The three lines show Y = 10, 1, 0.1 X. There is only one galaxy that haveṗout >ṗmax, which is similar to Figure 13. See Section 6.6 for discussion. locities of the warm phase of ∼ 10 3 km/s for all cases considered (roughly three to four times higher than our scaling relations in Figure 9 for an M82-like starburst). We conclude that these cooling wind models and numerical simulations are not a good match to our data.
Most recently, Lochhaas et al. (2021) examined the amount of momentum that can be carried by a warm outflow in the context of the contamination of the wind fluid (β > 1) and the resulting radiative losses. Recall that in the simple model described in Section 6.1, the amount of momentum flux that can be carried by the wind fluid isṗ wind = (αβ) 1/2ṗ * . Lochhaas et al. (2021) calculated the maximum allowed value of the product of αβ under the conditions of increasingly large β leading to increasing strong radiative cooling (α < 1). They find the maximum wind momentum flux to be: p max = 7.2 × 10 34 (R * ,0.3 ) 0.14 (αSFR 10 ) 0.86 dynes We compare our estimates of the outflow momentum flux withṗ max in Figure 16. We find that only one galaxy lies significantly above this relationship. In fact, this plot is very similar to our results in Figure 13, which simply compared the observed momentum flux to that provided by the starburst. The agreement means thatṗ max is generally similar tȯ p * . That is to say, actual outflows seem to carry the maximum possible momentum allowed for them without being quenched by radiative cooling, and this maximum flux is very similar to the amount input by the starburst (αβ ∼ 1).

CONCLUSION
We have reported here the results of our analysis of starburst-driven galactic outflows of warm ionized gas. This was based on data for 45 galaxies taken from the COS Legacy Archive Spectroscopic SurveY (CLASSY), augmented by five additional starbursts with COS data of similarly-high quality and under the same selection criteria. The properties of the outflows were based on fitting the UV resonance absorption-lines, and in particular, using five Si II multiplets and the Si IV 1393, 1402 doublet to derive the column density and covering fraction of these ions as a function of outflow velocity. CLOUDY models were used to derive total Si column densities and these were converted into H column densities using the metallicities derived from the nebular emission-lines.
The key parameters obtained from this analysis are the mean outflow velocity (V out ) and the Full Width at Half Maximum (FWHM out ) of the blue-shifted absorption-lines, the total Si and H column densities (N Si and N H ), and the outflow rates of Si (Ṁ Si,out ), mass (Ṁ out ), momentum (ṗ out ), and kinetic energy (Ė out ).
We then examine the scaling relationships between the outflow properties and those of the starburst and its host galaxy. The principle results are as follows: • Outflows were detected in roughly 90% of the sample. This implies that the outflows in starburst galaxies cover most of 4π steradian (they cannot be wellcollimated).
• The average value of the covering factor is 0.64, meaning that the effects of partial covering need to be explicitly determined to measure column densities from optically-thick absorption-lines. The values for the covering fractions are very similar for Si II and Si IV.
• While the values for V out are quite consistent among all the transitions we measure, there is significant scatter in the values for FWHM outİ n particular, we find a systematic trend for FWHM out to decrease as we move from the most to the least optically-thick Si II transitions. This implies that the highest column densities in the outflow are near the characteristic outflow velocity given by V out .
• We found highly significant correlations of both V out and FWHM out with star-formation rate (SFR), galaxy stellar mass (M ), and galaxy circular velocity (V cir ). The best-fit relationship is V out ∝ V cir 0.6 , and the ratio of V out /V cir being ∼ 6 and ∼ 2 for the lowest and highest mass galaxies, respectively.
• We found that the outflow rate of Si is (on-average) only about 20% of the rate at which Si is created and ejected by supernovae in the starburst. We conjecture that most of this "missing" Si is in the form of a hotter and more highly-ionized phase of the outflow than what we probe with these data.
• Assuming a normal dust-to-gas-phase-metals ratio, the observed Si column densities implied that there is very little dust extinction associated with the observed outflows (mean FUV optical depth of ∼ 0.3). Most of the observed extinction is produced by dust in the static ISM.
• The average total hydrogen column density (N H ) is ∼ 10 20.7 cm −2 , and that neutral hydrogen [N(H I)] only constitutes 0.1 to 1% of N H . The dominant ion for silicon is Si III, and 90 to 99% of the Si II arises in the ionized gas. Based on the derived N H (v), we find the column densities peak near V out , while the broad wings of outflow profiles have significantly smaller values of N H .
• We found a highly significant correlation between the mass outflow rate (Ṁ out ) and the SFR, but with a shallow slope (Ṁ out ∝ SFR 0.4 ). Hence the so-called massloading factor of the outflow isṀ out /SFR ∝ SFR −0.6 . We also found that the mass-loading factor is a steep inverse function of both V out and V cir (slope ∼ −1.6), with mass-loading factors of ∼ 10 in the lowest mass galaxies. Together with the third result above, this supports the idea of a mass-dependent impact of outflows on the evolution of galaxies.
• We found strong correlations between the rates at which the outflows carry momentum and kinetic energy (ṗ out andĖ out ) and the rates at which the starbursts supply momentum and kinetic energy (ṗ andĖ ). The median values areṗ out ∼ 30%ṗ andĖ out ∼ 5%Ė .
We then compared these results to various theoretical models of galactic winds driven by starbursts. We began with a description of the most widely-used model of galactic winds due to Chevalier & Clegg (1985). In this model the stellar ejecta in the starburst create very hot gas which flows out to create a supersonic wind. This fast-moving and tenuous "wind-fluid" interacts with ambient gas which it accelerates to create the warm ionized outflows we observe through the transfer of momentum. We then examined some specific models ranging from very simple analytic ones to state-ofthe-art hydrodyamical simulations.
• Following Heckman et al. (2015), we evaluated the ratio of the outward force on ambient gas clouds to the inward force of gravity. We found this ratio to have a median value of ∼ 30, with only 2 cases having a ratio < 10. Thus, this sample is in the regime of "strong outflows" in which gravitational forces are secondary.
• We considered a simple model of a wind-blown bubble driven by the momentum supplied by the starburst. While this model can fit the relationship between V out and SFR, it fails in other respects, and would produce absorption-lines with V out >> FWHM out (while we find V out ∼ 0.5 FWHM out ).
• We compared our results to a new semi-analytic model of multi-phase galactic winds. We found agreement, but only for models in which the hot wind-fluid that accelerates the warm ionized gas we observe is uncontaminated (made up of pure stellar ejecta) and the clouds it interacts with are massive (10 5 to 10 6 M ).
• Recent high-resolution numerical simulations by Schneider et al. (2020) produced outflows with some similarities to the data, but whose warm ionized phase was significantly too fast and carried too little mass, momentum, and kinetic energy compared to the data. In contrast, the simulations by Kim et al. (2020) predict the slope of outflow velocity vs. SFR/Area ∼ 0.2, which is similar to our observations. But their normalization of outflow velocity is small compared to the data.
• Finally, we compared the data to a family of models in which the warm ionized gas is not ambient material, but instead forms directly from the wind fluid (via radiative cooling). These models predicted outflow velocities significantly larger than we observed.
In a future paper, we will analyze FUSE spectra of starbursts using the same methodology as in this paper. This will allow us to extend the range of parameter space we can explore to lower values of both SFR/M * and SFR/Area, which will allow us to disentangle the dependence of outflow properties on SFR vs. M * vs. size. We will also analyze the current CLASSY spectra, searching for the possible presence of absorption lines arising from the fine structure levels in the Si II transitions. These could allow us to measure a mean density of the absorbing gas in the outflows for the first time.
The CLASSY team is grateful for the support for this program, HST-GO-15840, that was provided by NASA through a grant from the Space Telescope Science Institute, which is operated by the Associations of Universities for Research in Astronomy, Incorporated, under NASA contract NAS5-26555. BLJ thanks support from the European Space Agency (ESA). The CLASSY collaboration extends special gratitude to the Lorentz Center for useful discussions during the "Characterizing Galaxies with Spectroscopy with a view for JWST" 2017 workshop that led to the formation of the CLASSY collaboration and survey. This research has made use of the HSLA database, developed and maintained at STScI, Baltimore, USA.
* These objects are marked as "No Outflow" because less than half of their observed absorption troughs show blue-shifted outflow component. Note that the rest (< half) of their absorption troughs may still show outflow features. See discussion in Section 4.2.