# CLASSY III: 111 Based on observations made with the NASA/ESA Hubble Space
Telescope, obtained from the Data Archive at the Space Telescope Science
Institute, which is operated by the Association of Universities for Research
in Astronomy, Inc., under NASA contract NAS 5-26555. The Properties of
Starburst-Driven Warm Ionized Outflows

Xinfeng Xu Center for Astrophysical Sciences, Department of Physics &
Astronomy, Johns Hopkins University, Baltimore, MD 21218, USA  Timothy Heckman
Center for Astrophysical Sciences, Department of Physics & Astronomy, Johns
Hopkins University, Baltimore, MD 21218, USA  Alaina Henry Center for
Astrophysical Sciences, Department of Physics & Astronomy, Johns Hopkins
University, Baltimore, MD 21218, USA  Space Telescope Science Institute, 3700
San Martin Drive, Baltimore, MD 21218, USA  Danielle A. Berg Department of
Astronomy, The University of Texas at Austin, 2515 Speedway, Stop C1400,
Austin, TX 78712, USA  John Chisholm Department of Astronomy, The University
of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712, USA  Bethan
L. James AURA for ESA, Space Telescope Science Institute, 3700 San Martin
Drive, Baltimore, MD 21218, USA  Crystal L. Martin Department of Physics,
University of California, Santa Barbara, Santa Barbara, CA 93106, USA  Daniel
P. Stark Steward Observatory, The University of Arizona, 933 N Cherry Ave,
Tucson, AZ, 85721, USA  Alessandra Aloisi Space Telescope Science Institute,
3700 San Martin Drive, Baltimore, MD 21218, USA  Ricardo O. Amorín Instituto
de Investigación Multidisciplinar en Ciencia y Tecnología, Universidad de La
Serena, Raul Bitrán 1305, La Serena 2204000, Chile  Departamento de
Astronomía, Universidad de La Serena, Av. Juan Cisternas 1200 Norte, La Serena
1720236, Chile  Karla Z. Arellano-Córdova Department of Astronomy, The
University of Texas at Austin, 2515 Speedway, Stop C1400, Austin, TX 78712,
USA  Rongmon Bordoloi Department of Physics, North Carolina State University,
421 Riddick Hall, Raleigh, NC 27695-8202, USA  Stéphane Charlot Sorbonne
Université, CNRS, UMR7095, Institut d’Astrophysique de Paris, F-75014, Paris,
France  Zuyi Chen Steward Observatory, The University of Arizona, 933 N Cherry
Ave, Tucson, AZ, 85721, USA  Matthew Hayes Stockholm University, Department of
Astronomy and Oskar Klein Centre for Cosmoparticle Physics, AlbaNova
University Centre, SE-10691, Stockholm, Sweden  Matilde Mingozzi Space
Telescope Science Institute, 3700 San Martin Drive, Baltimore, MD 21218, USA
Yuma Sugahara Institute for Cosmic Ray Research, The University of Tokyo,
Kashiwa-no-ha, Kashiwa 277-8582, Japan  National Astronomical Observatory of
Japan, 2-21-1 Osawa, Mitaka, Tokyo 181-8588, Japan  Waseda Research Institute
for Science and Engineering, Faculty of Science and Engineering, Waseda
University, 3-4-1, Okubo, Shinjuku, Tokyo 169-8555, Japan  Lisa J. Kewley
Research School of Astronomy and Astrophysics, Australian National University,
Cotter Road, Weston Creek, ACT 2611, Australia; ARC Centre of Excellence for
All Sky Astrophysics in 3 Dimensions (ASTRO 3D), Canberra, ACT 2611, Australia
Masami Ouchi National Astronomical Observatory of Japan, 2-21-1 Osawa, Mitaka,
Tokyo 181-8588, Japan  Institute for Cosmic Ray Research, The University of
Tokyo, Kashiwa-no-ha, Kashiwa 277-8582, Japan  Kavli Institute for the Physics
and Mathematics of the Universe (WPI), University of Tokyo, Kashiwa, Chiba
277-8583, Japan  Claudia Scarlata Minnesota Institute for Astrophysics,
University of Minnesota, 116 Church Street SE, Minneapolis, MN 55455, USA
Charles C. Steidel Cahill Center for Astronomy and Astrophysics, California
Institute of Technology, MC249-17, Pasadena, CA 91125, USA Xinfeng Xu
<EMAIL_ADDRESS>
###### Abstract

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.

Galactic Winds (572), Galaxy evolution (1052), Galaxy kinematics and
dynamics(602), Starburst galaxies (1570), Ultraviolet astronomy (1736), Galaxy
spectroscopy (2171)

††journal: AASJournal ApJ††facilities: HST (COS), LBT (MODS), APO (SDSS), KECK
(ESI), VLT (MUSE, VIMOS)††software:  astropy (The Astropy Collaboration 2013,
2018) BEAGLE (Chevallard & Charlot 2006), CalCOS (STScI), dustmaps (Green
2018), jupyter (Kluyver 2016), MPFIT (Markwardt 2009), Photutils (Bradley
2021), python, pysynphot (STScI Development Team)

##  1 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 build-up 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. 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. 2.

Why does this efficiency reach its peak value over a relatively narrow range
in dark matter halo masses of $\sim 10^{12}$ M⊙, and why is this value roughly
constant with $z$?

  3. 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. 4.

Why is there such small scatter in the correlation between M⋆ and the galaxy’s
chemical composition (metallicity) at a given epoch?

  5. 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. 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 “sub-grid 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; 2020).

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; 2017;
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 (Berg et al.
2022; 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 signal-to-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 before, 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.

We adopt a cosmology with H0 = 69.6 km s-1 Mpc-1, $\Omega_{m}$ = 0.286, and
$\Omega_{\Lambda}$ = 0.714, and we use Ned Wright’s Javascript Cosmology
Calculator website (Wright 2006).

##  2 Observations

CLASSY is a Hubble Space Telescope (HST) treasury program (GO: 15840, PI:
Berg), which provides the first high-resolution 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⋆ $\sim 6-10\
M_{\odot}$), SFR ($\sim 0.01-100$ M⊙ yr-1), metallicity (12+log(O/H)$\sim
7-9$), and electron density ($n_{e}\sim 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.

##  3 Data Reduction

After the observations, all data were reduced locally using the COS data-
reduction package CalCOS
v.3.3.10222https://github.com/spacetelescope/calcos/releases. 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$ $\sim$ 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$ $\sim$ 6000 – 10000 from the blue
to red end) to gain higher S/N.

###  3.1 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.

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.

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 $\chi^{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 $\lambda$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.

###  3.2 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.

##  4 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 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.

###  4.1 Justifications of Methodology

For each galaxy in our sample, the final reduced, co-added, and starlight
subtracted spectra cover $\sim$ 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 $\lambda$1302, C
ii $\lambda$1334, and Si ii multiplet ($\lambda$1190, 1193, 1260, 1304, and
1526), and from higher-ionization transitions, e.g., Si iii $\lambda$1206, and
Si iv $\lambda$$\lambda$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 line-spread 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 absorption-line 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 2011). Second, it provides an excellent fit
to the data (as we will show).

Table 1: Atomic Data for Ions Measured in Outflows(∗) Ions | Vac. Wave. |  flk |  Akl |  El – Ek
---|---|---|---|---
(1) | (2) | (3) | (4) | (5)
O i | 1302.17 |  5.20 $\times$ 10-2 |  3.41 $\times$ 108 | 0.0 - 9.52
C ii | 1334.53 |  1.29 $\times$ 10-1 |  2.41 $\times$ 108 | 0.0 - 9.29
Si ii | 1190.42 |  2.77 $\times$ 10-1 |  6.53 $\times$ 108 | 0.0 - 10.41
Si ii | 1193.29 |  5.75 $\times$ 10-1 |  2.69 $\times$ 109 | 0.0 - 10.39
Si ii | 1260.42 |  1.22 $\times$ 10-1 |  2.57 $\times$ 109 | 0.0 - 9.84
Si ii | 1304.37 |  9.28 $\times$ 10-2 |  3.64 $\times$ 108 | 0.0 - 9.50
Si ii | 1526.71 |  1.33 $\times$ 10-1 |  3.81 $\times$ 108 | 0.0 - 8.12
Si iii | 1206.51 | 1.67 |  2.55 $\times$ 109 | 0.0 - 10.27
Si iv | 1393.76 |  5.13$\times$ 10-1 |  8.80 $\times$ 108 | 0.0 - 8.90
Si iv | 1402.77 |  2.55$\times$ 10-1 |  8.63 $\times$ 108 | 0.0 - 8.84
Note. – (*). Data obtained from National Institute of Standards and Technology
(NIST) atomic database (Kramida et al. 2018). (2). Vacuum wavelengths in units
of  Å. (3). Oscillator strengths. (4). Einstein A coefficients in units of
s-1. (5). Energies from lower to higher levels in units of eV.

###  4.2 Double-Gaussian Fits of the Absorption Troughs

The theoretical LSF provided in HST/COS
website333https://www.stsci.edu/hst/instrumentation/cos/performance/spectral-
resolution is good for a point source, but CLASSY galaxies are largely
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. 1.

First of all, we would like to generate a LSF for each galaxy based on the COS
LSF for a point source (LSF0) 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
(FWHMuv, see Table 2 in Paper I). An example is shown in Figure 1). We assume
the galaxy has Gaussian profile with FWHMuv (Guv, in blue), which is
sufficient for our analyses. Then we convolve LSF0 (in red) with Guv, which
results in an approximate non-point source LSF (LSFuv, in black) for a given
galaxy. This LSFuv is then used in the double-Gaussian fittings below to
properly account for the LSF of each galaxy.

  2. 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
LSFuv 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. 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:

$F=\frac{(\chi^{2}_{1}-\chi^{2}_{2})/(p_{2}-p_{1})}{\chi^{2}_{2}/(n-p_{2})}$
(1)

where $\chi^{2}_{1}$ and $\chi^{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 $\alpha$ =
0.05 and degree of freedom of ($p_{2}$ – $p_{1}$, $n$ – $p_{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.

Figure 3: Distributions for outflow velocity (Left) and the FWHM of the
outflow component (Right) for galaxies in our combined sample [CLASSY +
Heckman et al. (2015), see Section 3]. The method for estimating these two
values are discussed in Section 4.2. For FWHM${}_{\text{out}}$, we have
already subtracted the contributions from the HST/COS line-spread-functions
(LSF).

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 LSFuv are usually $\ll$ 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${}_{\text{out}}$):

$\text{FWHM}_{\text{out}}=\sqrt{\text{FWHM}_{\text{all}}^{2}-\text{FWHM}_{\text{LSF}}^{2}}$
(2)

where $\text{FWHM}_{\text{all}}$ is the fitted FWHM for a certain trough and
$\text{FWHM}_{\text{LSF}}$ is the FWHM of LSFuv for that galaxy.

For each galaxy, we adopt this method to fit the absorption lines from O i
$\lambda$1302, C ii $\lambda$1334, Si ii multiplet ($\lambda$1190, 1193, 1260,
1304, and 1526), Si iii $\lambda$1206, and Si iv $\lambda$$\lambda$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$\alpha$), 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${}_{\text{out}}$) is derived from the median value of FWHM from all
troughs that pass the F-test. The corresponding errors of $V_{\text{out}}$ (or
FWHM${}_{\text{out}}$) are estimated from the standard deviations of
$V_{\text{out}}$ (or FWHM${}_{\text{out}}$) from all passed lines. Note that
for these “no outflow” galaxies, they could still host very low velocity
outflows ($V_{\text{out}}$$\ll$ FWHM${}_{\text{ISM}}$). However, since we
cannot disentangle them from the static ISM component, their $V_{\text{out}}$
and FWHM${}_{\text{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 $\lambda$1304
(which is the most optically thin Si ii line we measure) to be narrower than
that of Si ii $\lambda$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 ($\sim$ 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_{\text{out}}$ and FWHM${}_{\text{out}}$ are shown in Figure 3, and the
values are presented Table 5.

Figure 4: Comparisons of outflow velocity centers ($V_{\text{out}}$) and
widths (FWHM${}_{\text{out}}$) from Si ii $\lambda$1260 and Si iv
$\lambda$1393 troughs (top) and for Si ii $\lambda$1260 and Si ii
$\lambda$1304 (bottom) in our combined sample. While all three lines show
consistent $V_{\text{out}}$, their FWHM${}_{\text{out}}$ values have
substantial scatter. See Section 4.2 for details. The black lines show the 1:1
relationship, and the orange lines show the best linear fits while the fitted
parameters are shown in the left-top corners. The results of Pearson
correlation coefficients (PCC) are shown at the bottom-right corners.

Figure 5: Examples of checking the amount of infilling by scattered resonance
emission based on the properties of the associated fluorescent emission-lines.
For double-Gaussian fits, the blue, green, and red lines are for the outflow,
static ISM, and summation of both components, respectively. The emission from
fluorescent emission-lines are fitted by one Gaussian profile shown as the
orange lines. For three-Gaussian fits, the orange dashed lines centered at $v$
= 0 km s-1 represent the amount of infilling, which are scaled and shifted
from the fluorescent lines. Left: A case when the infilling is negligible.
Right: A case when the infilling is obvious at the bottom of static ISM
component. This infilling will only affect the fit of the static ISM
component, but not the outflow component. See details in Section 4.3.

###  4.3 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*
$\lambda$1265 for Si ii $\lambda$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 $\lambda$1260 would
cover the same velocity range as Si ii* $\lambda$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 $\lambda$$\lambda$1393, 1402, but there are no associated fluorescent
transitions for them. The steps are as follows:

  1. 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. 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 $\lambda$1260, this
constant equals $A_{1260}/A_{1265}\times H_{1260}\simeq 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. 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).

Figure 6: An example of the solutions of partial-covering equations and total
hydrogen column densities ($N_{\mathrm{\scriptstyle H}}$) 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 $\lambda$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 $\lambda$1304 line is
significantly narrower than the other three lines Top-Right: Best fitting
covering fraction (C${}_{\text{f}}$) for Si ii from solving Equation (4).
Bottom-Left: Best fitting optical depth over velocity ($\tau$(v)) in logarithm
for each Si ii lines. We only show Si ii $\lambda$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
C${}_{\text{f}}$(v)$\times\tau$(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.

###  4.4 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_{\text{cir}}$) and r50. Following
Heckman et al. (2015), we calculate $V_{\text{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_{\text{cir}}$ $=\sqrt{2}$ $(V_{rot}^{2}+2\sigma^{2})^{1/2}$, where
$V_{rot}$ and $\sigma$ are the rotation speed and mean value of the line-of-
sight velocity dispersion, respectively. Finally, we adopt the best-fit
relation from Simons et al. (2015) as log $V_{\text{cir}}$ = 0.29 M⋆ – 0.78,
with root-mean-square (RMS) residuals of 0.1 dex RMS in $V_{\text{cir}}$.

For each galaxy, we measure the half-light radius (r50) 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 r50 values from COS for
galaxies with r50 smaller than 0.4″. These galaxies are compact enough, so the
vignetting effects are small. For other galaxies that have COS r50 $>$ 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 r50 $>$ 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 instead 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_{\text{cir}}$ and r50 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.

##  5 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${}_{\text{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_{\mathrm{\scriptstyle 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.

###  5.1 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${}_{\text{f}}$) $<$ 1.0,
the column density ($N_{\mathrm{\scriptstyle 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_{\mathrm{\scriptstyle ion}}$ and
C${}_{\text{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; 2017). 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${}_{\text{f}}$ and the optical depth ($\tau$) 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
$\lambda$$\lambda$1393, 1402), we have:

$\displaystyle I_{R}(v)$
$\displaystyle=1-C_{\text{f}}(v)+C_{\text{f}}(v)\times e^{-\tau(v)}$ (3)
$\displaystyle I_{B}(v)$
$\displaystyle=1-C_{\text{f}}(v)+C_{\text{f}}(v)\times e^{-w\cdot\tau(v)}$

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 $\tau(v)$ is the optical depth of the red line (Si iv $\lambda$1402) that
is velocity dependent. The weight $w$ in the bottom equation equals
$f_{B}\lambda_{B}$/$f_{R}\lambda_{R}$, where $f$ is the oscillator strength
and $\lambda$ is the wavelength of the line (see Table 1). For common doublet
lines, such as N v $\lambda$$\lambda$1238, 1242, C iv $\lambda$$\lambda$1548,
1550, and Si iv $\lambda$$\lambda$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:

$\displaystyle I_{k}(v)$
$\displaystyle=1-C_{\text{f}}(v)+C_{\text{f}}(v)\times e^{-w_{k}\cdot\tau(v)}$
(4)

where for Si ii multiplet, $k$ stands for Si ii $\lambda$1190, 1193, 1260,
1304, and 1526. We define $w_{k}$ =
$f_{k}\lambda_{k}$/$f_{1304}\lambda_{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 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 $\lambda$$\lambda$ 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 $\tau(v)$ and
C${}_{\text{f}}(v)$ from Equations (3) and (4) for Si iv and Si ii,
respectively. An example for the fitted Si ii troughs and the resulting
C${}_{\text{f}}(v)$, $\tau(v)$, and $N_{\mathrm{\scriptstyle ion}}$ are shown
in Figure 6. For each galaxy, we also list the outflow mean covering fraction
($\overline{C_{\text{f}}}$) in Table 5. This is calculated from the average of
C${}_{\text{f}}(v)$ in the range of [$V_{\text{out}}$ \-
$\sigma_{\text{Vout}}$,$V_{\text{out}}$ \+ $\sigma_{\text{Vout}}$] for both
the Si ii and Si iv multiplet/doublet lines. This choice is because most of
the outflowing material (i.e., $N_{\mathrm{\scriptstyle ion}}$) has velocities
around the center of the troughs. Therefore, $\overline{C_{\text{f}}}$
represents the covering fraction for the bulk of outflowing material we
observe from Si ii and Si iv.

Finally, we calculate $N_{\mathrm{\scriptstyle ion}}$ for Si ii and Si iv as
below (e.g., Savage & Sembach 1991; Edmonds et al. 2011).

$\displaystyle\bar{N}_{ion}(v)$ $\displaystyle=\frac{3.8\times\
10^{14}}{f\cdot\lambda}\cdot C_{\text{f}}(v)\cdot\tau(v)$ (5) $\displaystyle
N_{ion}$ $\displaystyle=\int{\bar{N}_{ion}(v)}dv$

where $\bar{N}_{ion}(v)$ is the average column density per velocity over the
aperture at $v$, and $N_{\mathrm{\scriptstyle 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 $\lambda$ 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 $\lambda$ 1206 is close to Ly$\alpha$, where the damped Ly$\alpha$
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 $\lambda$ 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 $\lambda$ 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
$\lambda$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_{\text{out}}$ and FWHM${}_{\text{out}}$ values from Si
ii $\lambda$1260 and 1304 for all our galaxies in the bottom two panels of
Figure 4. We find troughs from Si ii $\lambda$1260 indeed commonly have larger
FWHM${}_{\text{out}}$ than that of Si ii $\lambda$1304, while their
$V_{\text{out}}$ are quite consistent. These results imply that the
FWHM${}_{\text{out}}$ parameter is quite sensitive to the distribution of
column densities of the outflows (unlike $V_{\text{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.

Figure 7: Zoom-in on the spectral regions around S ii
$\lambda$$\lambda$$\lambda$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 $\lambda$1260\. Therefore, we exclude it when
calculating the column density of S ii (see Section 5.2).

###  5.2 Dust Depletion of Gas-Phase Metals

To estimate the total hydrogen column density ($N_{\mathrm{\scriptstyle 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); Jenkins (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
$\lambda$$\lambda$$\lambda$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 SNR2
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$\sigma$, $med$ \+ 3$\sigma$], where
$med$ and $\sigma$ 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.

Table 2: Column Densities for the Stacked Spectra Labels |  # of Galaxy(1) |  N(Si ii)${}_{\text{obs}}$ |  N(S ii)${}_{\text{obs}}$ |  N(S ii)${}_{\text{pred}}$(2)
---|---|---|---|---
S ii stack | 31 |  774${}^{+12}_{-12}$ |  168${}^{+35}_{-35}$ |  203${}^{+5}_{-5}$
|  |  N(Si iv) |  N(S iv)${}_{\text{obs}}$ |  N(S iv)${}_{\text{pred}}$(2)
S iv stack | 24 |  500${}^{+3}_{-3}$ |  193${}^{+39}_{-39}$ |  175${}^{+2}_{-2}$
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)).

We follow the same methodology in Section 5.1 to measure the column densities
of S ii [i.e., N(S ii)${}_{\text{obs}}$]. For S ii $\lambda$1259, since it is
too close to Si ii $\lambda$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:

$\displaystyle\text{ N(S II})_{\text{pred}}$ $\displaystyle=\mkern
1.5mu\overline{\mkern-1.5mu\text{N}\mkern-1.5mu}\mkern 1.5mu\text{(Si
II})\cdot(\text{S}/\text{Si})\cdot\frac{\text{ICF(S II})}{\text{ICF(Si II})}$
(6)

where S/Si is the abundance ratio, and we assume solar abundance; $\mkern
1.5mu\overline{\mkern-1.5mu\text{N}\mkern-1.5mu}\mkern 1.5mu$(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) 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) $\sim$ 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)${}_{\text{obs}}$ and N(S
ii)${}_{\text{pred}}$ from our stacked spectra are consistent within
1$\sigma$. Similarly, we conduct another stack which includes all galaxies
covering S iv $\lambda$1062 region, and find N(S iv)${}_{\text{obs}}$ is also
consistent with N(S iv)${}_{\text{pred}}$ within 1$\sigma$ (see the second
part of Table 2). In this case, N(S iv)${}_{\text{pred}}$ is scaled from
$\mkern 1.5mu\overline{\mkern-1.5mu\text{N}\mkern-1.5mu}\mkern 1.5mu$(Si iv)
given a solar S/Si ratio and ICF(S iv)/ICF(Si iv) $\sim$ 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. Qualitatively-similar 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_{\mathrm{\scriptstyle H}}$ of outflows from the measured N(Si ii) and N(Si
iv), which is discussed next.

###  5.3 Column Densities from CLOUDY Models

####  5.3.1 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_{\mathrm{\scriptstyle Si}}$) and total
hydrogen column density ($N_{\mathrm{\scriptstyle 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_{\mathrm{\scriptstyle e}}$) = 10 cm-3 for
typical starburst 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_{\mathrm{\scriptstyle Si}}$ to
$N_{\mathrm{\scriptstyle H}}$ assumes that the outflow has the same
metallicity as measured from the nebular emission lines.444In principle, the
metallicity of the outflow could be larger than in the HII regions, if the
former is strongly contaminated by metals ejected by the supernovae (Chisholm
et al. 2018; Hogarth et al. 2020). We will show in Section 5.4.3 that this is
unlikely. We have tested models for $n_{\mathrm{\scriptstyle e}}$ = 100 cm-3,
and find that the resulting $N_{\mathrm{\scriptstyle Si}}$ and
$N_{\mathrm{\scriptstyle H}}$ only have minor changes ($<$ 3 – 5%). In our
grid models, the two varied parameters are: 1) the logarithm of ionization
parameter, log($U_{\mathrm{\scriptstyle H}}$), in the range between -4.0 and
0.0 with a step size of 0.05 dex, and 2) the logarithm of
$N_{\mathrm{\scriptstyle 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_{\mathrm{\scriptstyle H}}$($v$) and $N_{\mathrm{\scriptstyle 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 $\chi^{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_{\mathrm{\scriptstyle H}}$($v$)
over all velocity bins to get $N_{\mathrm{\scriptstyle H}}$. Hereafter, we use
$N_{\mathrm{\scriptstyle ion}}$($v$) to represent the column density per
velocity at $v$ [in unit of cm-2/(km s-1)], while $N_{\mathrm{\scriptstyle
ion}}$ stands for the integrated column density over all velocity bins in the
outflow trough [in unit of cm-2]. The best-fit $N_{\mathrm{\scriptstyle H}}$
values are listed in Table 5, and its distribution is shown in Figure 8, where
we find that the mean $N_{\mathrm{\scriptstyle H}}$ for the observed outflows
is 1020.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_{\mathrm{\scriptstyle Si}}$ = N(Si ii) +
N(Si iii) + N(Si iv).

Furthermore, the mean value of N(Si iii) is 1015.53 cm-2. For an average
FWHMout $\sim$ 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:

$\tau\ =\frac{N_{\text{Si III}}\cdot f_{\text{Si III}}\cdot\lambda_{\text{Si
III}}}{3.8\times\ 10^{14}\cdot FWHM_{out}}$ (7)

where we have $f_{\text{Si III}}$ = 1.67 and $\lambda_{\text{Si III}}$ =
1206.51 Å. This leads to $\tau$(Si iii) $\sim$ 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 outflows 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).

####  5.3.2 The Neutral Phase of the Outflow

For our galaxies, the mean N(H i) from CLOUDY models is 1018.28 cm-2, which is
only $\sim$ 0.4% of the mean $N_{\mathrm{\scriptstyle H}}$ value (1020.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 $\lambda$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.

###  5.4 Scaling Relationships

In this subsection, we present various empirical scaling relationships between
the outflow and galaxy properties. These correlations for low-redshift
galaxies (z $\lesssim$ 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 low- and 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; 2017; 2018). 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.

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., $V_{\text{out}}$ $\pm$ $\sigma_{\text{Vout}}$ (see Section 5.1). Right:
Distributions for the total column density of outflows
($N_{\mathrm{\scriptstyle H}}$, 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.

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.

Figure 9: The log of the outflow velocity ($V_{\text{out}}$) and Full-Width-
Half-Maximum (FWHM${}_{\text{out}}$) plotted as a function of the basic
properties of the star-forming galaxies. From top-left to bottom-right are
log($V_{\text{out}}$) vs. circular velocity, star formation rate, and
log(FWHM${}_{\text{out}}$) 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
$V_{\text{out}}$$=V_{cir}$ and $V_{\text{out}}$$=10V_{cir}$ (and the same for
FWHM${}_{\text{out}}$). Galaxies with non-detections of outflows are shown as
hollow circles. We set their $V_{\text{out}}$ = 40 km s-1 and
log(FWHM${}_{\text{out}}$) = 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)].

Figure 10: Comparisons of normalized outflow velocity
($V_{\text{out}}$/$V_{\text{cir}}$) 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 $V_{\text{out}}$/$V_{\text{cir}}$ in our sample,
the correlations in these two panels are weak (L) to insignificant (R). See
discussion in Section 5.4.1.

####  5.4.1 Outflow Velocity and FWHM

We begin by examining the correlations of $V_{\text{out}}$ and
FWHM${}_{\text{out}}$ (Section 4.2) with the principal properties of the
galaxy, including the circular velocity ($V_{\text{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_{\text{out}}$ =
40 km s-1 and FWHM${}_{\text{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 $\tau$ 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_{\text{out}}$ (or FWHM${}_{\text{out}}$) with M⋆ (and
hence with $V_{\text{cir}}$). It is noteworthy that the slopes of the
relationships between $V_{\text{out}}$ and FWHM${}_{\text{out}}$ with
$V_{\text{cir}}$ are sub-linear, meaning that their ratios decrease with
increasing $V_{\text{cir}}$. We will discuss this further in Section 6 below.
In the right column of Figure 9, we see strong positive correlations between
$V_{\text{out}}$ (or FWHM${}_{\text{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_{\text{out}}$ with
SFR than with M⋆.

We have also noticed that the correlations with FWHM${}_{\text{out}}$ are
slightly stronger than the ones with $V_{\text{out}}$ in Figure 9. As we
explained in Section 5.1, unlike $V_{\text{out}}$, FWHM${}_{\text{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${}_{\text{out}}$ somewhat differently than $V_{\text{out}}$.

In Figure 10, we also test the correlations between the normalized outflow
velocity ($V_{\text{out}}$/$V_{\text{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_{\text{out}}$ to SFR, M⋆, and $V_{\text{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_{\text{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_{\text{out}}$, we only list the log-log 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_{\text{cir}}$, and M⋆.

Figure 11: Correlations related to the silicon mass flow rate ($\dot{M}_{\text{Si,out}}$) and total silicon column density ($N_{\mathrm{\scriptstyle Si}}$). Labels and captions are the same as Figure 9. Left: We compare $\dot{M}_{\text{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 $N_{\mathrm{\scriptstyle Si}}$ 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 $\sim$ 16% of the total dust measured. The right y-axis represents the dust optical depth at 1200 Å assuming Reddy et al. (2015) extinction law. See discussion in Section 5.4.3.  Table 3: Comparisons of $V_{\text{out}}$ Correlations with Published Star-forming/Starburst Galaxy Samples References | Redshift |  Lines(1) |  Def. of $V_{\text{out}}$(2) | SFR | $V_{\text{cir}}$ |  M⋆
---|---|---|---|---|---|---
This paper | 0.002 – 0.23 |  Far-UV lines(3) | Double Gaussian |  0.22$\pm$0.08 |  0.66$\pm$0.21 |  0.19$\pm$0.06
Martin (2005) | 0.04 – 0.16 |  Na i $\lambda$$\lambda$5890, 5896 | Double Gaussian |  0.35$\pm$0.06 | … | …
Rupke et al. (2005) | 0.01 – 0.50 |  Na i $\lambda$$\lambda$5890, 5896 | Max velocity |  0.21$\pm$0.05 | … | …
Chisholm et al. (2015)(4) | 0.0007 – 0.26 |  Si ii 1190, 1193, 1260, 1304 | Single Gaussian |  0.22$\pm$0.04 |  0.87$\pm$0.17 |  0.20$\pm$0.05
Chisholm et al. (2016a)(4) | 0.0007 – 0.26 |  Far-UV lines(3) | Max velocity |  0.12$\pm$0.02 | … |  0.15$\pm$0.02
Heckman & Borthakur (2016) | 0.4 – 0.7 |  UV lines(3) | Max velocity |  0.32$\pm$0.02 |  1.16$\pm$0.37 |  0.34$\pm$0.11
Note. – 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 $V_{\text{out}}$ and SFR, $V_{\text{cir}}$, and M⋆ are shown in
the 5th, 6th, and 7th columns, respectively. (1). The lines adopted to measure
the outflow velocity ($V_{\text{out}}$). (2). Definitions for
$V_{\text{out}}$. In this paper, we derive $V_{\text{out}}$ 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 $v_{90}$) of the troughs as $V_{\text{out}}$. (3). From HST/COS
spectra, the major rest-frame Far-UV lines adopted to estimate
$V_{\text{out}}$ 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 $\lambda$977, C ii
$\lambda$1036, and N ii $\lambda$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).

Figure 12: Correlations related to total mass flow rate
($\dot{M}_{\text{out}}$) and mass loading factor ($\dot{M}_{\text{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 $\dot{M}_{\text{out}}$/SFR scales with
the galaxy circular velocity and with the outflow velocity. In both cases, the
ratio increases rapidly at lower velocities (slopes $\sim$ –1.6 to –1.7). See
further discussion in Section 5.4.4.

####  5.4.2 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 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
($\dot{M}_{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 to
$\dot{M}_{out}=M_{out}V_{out}/R_{out}$, where $V_{\text{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
$\Omega$ we have:

$\displaystyle\dot{M}_{out}\simeq\Omega N_{\text{H}}\mu
m_{p}R_{\text{shell}}V_{\text{shell}}$ (8)

where $\mu$ is the mean mass per proton ($\sim 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_{\text{shell}}$. We will return to consideration of shells in the
context of a wind-blown bubble model in Section 6 below.

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_{\text{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 that

$\displaystyle\dot{M}_{out}\simeq\Omega N_{\text{H}}\mu
m_{p}R_{\text{0}}V_{\text{out}}$ (9)

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):

$\displaystyle\dot{M}_{out}\simeq\Omega N_{\text{H}}\mu
m_{p}R_{0}V_{max}/ln(R_{max}/R_{0})$ (10)

where $R_{max}$ is the maximum radius which the outflow reaches.

The key point here is that for reasonable choices for the relevant parameters,
the values of $\dot{M}_{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}$ $\sim$ $V_{\text{out}}$\+ 0.5 FWHM${}_{\text{out}}$, then Figure 5
implies $V_{max}$ $\sim$ 2 $V_{\text{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}$.

####  5.4.3 Metal Mass Outflow Rates and Dust Extinction

With these considerations in mind, we calculate the metal mass outflow rates
of silicon ($\dot{M}_{\text{Si,out}}$) for galaxies in our sample. One
advantage of $\dot{M}_{\text{Si,out}}$ over the total (i.e., hydrogen) mass
outflow rate ($\dot{M}_{\text{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)] for silicon:

$\displaystyle\dot{M}_{\text{Si,out}}\simeq\Omega\mu
m_{p}R_{\text{0}}\int\frac{dN_{\text{Si}}}{dv}\times vdv$ (11)

where $dN_{\text{Si}}/dv$ = N${}_{\text{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., $\Omega$,
$N_{\mathrm{\scriptstyle Si}}$, and $R_{0}$) required to calculate
$\dot{M}_{\text{Si,out}}$. For $\Omega$, as shown in Section 4.2, in the
CLASSY sample, $\sim$ 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 $\Omega$ =
4$\pi$.

The value for $N_{\mathrm{\scriptstyle 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$ $\times$
r50 (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 $\dot{M}_{\text{Si,out}}$ with the
rate at which silicon is created and injected by SNe in the starburst, i.e.,
$\dot{M}_{\text{Si,*}}$ (both in units of M⊙/yr). We approximate
$\dot{M}_{\text{Si,*}}$ = 1$\times$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 $\dot{M}_{\text{Si,out}}$ and
$\dot{M}_{\text{Si,*}}$ is only $\sim$ 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_{\mathrm{\scriptstyle Si}}$ and dust extinction are related to
the total column density of metals. The difference is that
$N_{\mathrm{\scriptstyle 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_{\mathrm{\scriptstyle
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_{\mathrm{\scriptstyle Si}}$, assuming a
standard Milky Way (MW) dust/metals ratio (Draine 2011) and assuming the solar
ratio of silicon to metals mass. This then implies that the total dust
extinction is related to $N_{\mathrm{\scriptstyle Si}}$ as (see Draine 2011):

$N_{\text{Si}}>E(B-V)_{\text{tot}}\times 1.9\times 10^{17}cm^{-2}$ (12)

In our sample, we then check the distribution of $N_{\mathrm{\scriptstyle
Si}}$/E(B-V)${}_{\text{tot}}$/(1.9 $\times$ 1017), and the resulting median is
$\sim$ 0.16. This suggests that the dust in our observed outflows only
represents $\sim$ 16% of E(B-V)${}_{\text{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
($\tau_{1200}$) at 1200 Å assuming the Reddy et al. (2015) extinction law and
Equation (12). This leads to a typical value of $<\tau_{1200}>$ $\sim$ 0.3 due
to dust in the outflow.

####  5.4.4 Total Mass Outflow Rates and Mass Loading Factors

Similarly to the calculations of $\dot{M}_{\text{Si,out}}$ [see Equation
(11)], we can derive the total mass outflow rate ($\dot{M}_{\text{out}}$) as:

$\displaystyle\dot{M}_{out}\simeq\Omega\mu
m_{p}R_{0}\int\frac{dN_{\text{H}}}{dv}\times vdv$ (13)

where $v$ is the outflow velcity, $dN_{\text{H}}/dv$ =
$N_{\mathrm{\scriptstyle 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 $\Omega$. In Figure 12, we present the correlations
related to $\dot{M}_{\text{out}}$.

In the first panel, we show a strong correlation between
$\dot{M}_{\text{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, $\dot{M}_{\text{out}}$ ranges from 0.3 –
40 M⊙/yr, and the best linear fit yields that $\dot{M}_{\text{out}}$ = 3.8
$\times$ SFR0.41 (a sub-linear slope).

In the next two panels, we show the correlations between the mass loading
factor (i.e., $\dot{M}_{\text{out}}$/SFR) and $V_{\text{cir}}$ (and
$V_{\text{out}}$). Both figures show strong inverse correlations, with the
correlation between $\dot{M}_{\text{out}}$/SFR and $V_{\text{cir}}$ being
stronger. The best fit slopes for both figures are $\sim-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 $\sim 10$ in the lowest mass
galaxies.

Figure 13: Correlations comparing (top left) the estimated momentum outflow
rate ($\dot{p}_{\text{out}}$) and (top right) the estimated kinetic energy
outflow rate ($\dot{E}_{\text{out}}$) to the rates of momentum and kinetic
energy supplied by the starburst. Labels and captions are the same as Figure
9. The three black lines refer to Y = 0.1, 1, 10 X in the left panel and Y =
0.01, 0.1, 1 X in the right panel. See discussion in Section 5.4.5. In the
last panel, we present comparisons to a simple analytic model of momentum-
driven outflows from Heckman et al. (2015) (see Section 6.2).

####  5.4.5 Momentum and Kinetic Energy Outflow Rates

Given the calculated $\dot{M}_{\text{out}}$, we can derive the momentum and
kinetic energy outflow rates as:

$\displaystyle\dot{p}_{out}$ $\displaystyle=\int\dot{M}_{out}(v)v=\Omega\mu
m_{p}R_{0}\int\frac{dN_{\text{H}}}{dv}\times v^{2}dv$ (14)
$\displaystyle\dot{E}_{out}$
$\displaystyle=\int\frac{1}{2}\dot{M}_{out}(v)v^{2}=\frac{1}{2}\Omega\mu
m_{p}R_{0}\int\frac{dN_{\text{H}}}{dv}\times v^{3}dv$

The integration over $v$ in $\dot{p}_{\text{out}}$ and $\dot{E}_{\text{out}}$
is necessary since both functions are strongly dependent on $v$. We can then
compare $\dot{p}_{\text{out}}$ and $\dot{E}_{\text{out}}$ with the momentum
and kinetic energy supplied from starburst, i.e., $\dot{p}_{\star}$ and
$\dot{E}_{\star}$. As discussed in Heckman et al. (2015), $\dot{p}_{\star}$ 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 to $\dot{p}_{\star}$ =
4.6$\times$1033 SFR dynes, where SFR is in units of M⊙/yr. Similarly,
$\dot{E}_{\star}$ = 4.3$\times$1041 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
$\dot{p}_{\text{out}}$ = 10 – 100% $\dot{p}_{\star}$, and
$\dot{E}_{\text{out}}$ = 1 – 20% $\dot{E}_{\star}$. These ranges are similar
to the ones reported in chisholm17, 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 $\dot{E}_{\text{out}}$$>$
$\dot{E}_{\star}$. However, there are three galaxies in which
$\dot{p}_{\text{out}}$$>$ $\dot{p}_{\star}$. 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.

Table 4: Fitted Parameters for Scaling Relationships Correlations |  r${}^{(1)}_{\text{k}}$ |  p${}^{(1)}_{\text{k}}$ |  r${}^{(2)}_{\text{PCC}}$ |  p${}^{(2)}_{\text{PCC}}$ |  Slope(3) |  Intercept(3)
---|---|---|---|---|---|---
FWHM${}_{\text{out}}$ vs. $V_{\text{cir}}$ | 0.55 | 8.88E-06 | 0.79 | 1.19E-07 |  0.52$\pm$0.22 |  1.62$\pm$0.38
FWHM${}_{\text{out}}$ vs. SFR | 0.62 | 7.15E-07 | 0.80 | 0.00E+00 |  0.19$\pm$0.08 |  2.38$\pm$0.09
$V_{\text{out}}$ vs. $V_{\text{cir}}$ | 0.53 | 2.15E-05 | 0.73 | 2.50E-06 |  0.65$\pm$0.21 |  1.14$\pm$0.37
$V_{\text{out}}$ vs. SFR | 0.49 | 8.70E-05 | 0.67 | 2.77E-05 |  0.22$\pm$0.08 |  2.12$\pm$0.09
$V_{\text{out}}$/$V_{\text{cir}}$ vs. sSFR | 0.30 | 1.64E-02 | 0.34 | 5.61E-02 |  0.16$\pm$0.16 |  1.84$\pm$1.28
$V_{\text{out}}$/$V_{\text{cir}}$ vs. SFR/Area | –0.18 | 1.44E-01 | –0.30 | 9.19E-02 |  –0.08$\pm$0.11 |  0.61$\pm$0.12
$\dot{M}_{\text{out}}$ vs. SFR | 0.29 | 3.29E-02 | 0.60 | 6.67E-04 |  0.38$\pm$0.10 |  0.64$\pm$0.11
$\dot{M}_{\text{out}}$/SFR vs. $V_{\text{out}}$ | –0.37 | 5.03E-03 | –0.59 | 9.04E-04 |  –1.73$\pm$0.43 |  4.07$\pm$0.95
$\dot{M}_{\text{out}}$/SFR vs. $V_{\text{cir}}$ | –0.52 | 9.16E-05 | –0.73 | 9.66E-06 |  –1.63$\pm$0.35 |  2.98$\pm$0.60
$\dot{p}_{\text{out}}$ vs. $\dot{p}_{\star}$ | 0.38 | 4.44E-03 | 0.71 | 1.97E-05 |  0.55$\pm$0.12 |  15.17$\pm$4.05
$\dot{E}_{\text{out}}$ vs. $\dot{E}_{\star}$ | 0.46 | 6.78E-04 | 0.75 | 5.13E-06 |  0.72$\pm$0.14 |  10.77$\pm$5.83
$\dot{M}_{\text{Si,out}}$/$\dot{M}_{\text{Si,*}}$ vs. $V_{\text{out}}$ | –0.18 | 1.75E-01 | –0.38 | 4.75E-02 |  –0.88$\pm$0.47 |  1.45$\pm$1.06
$\dot{M}_{\text{Si,out}}$/$\dot{M}_{\text{Si,*}}$ vs. $V_{\text{cir}}$ | –0.38 | 5.56E-03 | –0.57 | 1.95E-03 |  –0.97$\pm$0.39 |  1.16$\pm$0.68
N(Si) vs. E(B-V) | 0.19 | 1.65E-01 | 0.41 | 3.98E-02 |  –0.97$\pm$0.39 |  1.16$\pm$0.68
$\dot{M}_{\text{Si,out}}$ vs. $\dot{M}_{\text{Si,*}}$ | 0.47 | 5.82E-04 | 0.79 | 1.19E-06 |  0.58$\pm$0.11 |  –1.44$\pm$0.25
Note. – (1). Coefficients for Kendall $\tau$ test. (2). Coefficients for
Pearson Correlation (PCC). (3). Fitted slope and intercept assuming linear
correlations between the X and Y values. Note that except for “N(Si) vs.
E(B-V)”, all other correlations are fitted in log-log scale. See Figures 9 –
13.

##  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.

###  6.1 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
($\dot{M}_{\star}$ and $\dot{E}_{\star}$) 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:

$\displaystyle v_{wind}=(2\dot{E}_{wind}/\dot{M}_{wind})^{1/2}$ (15)

$\displaystyle\dot{M}_{wind}=\beta\dot{M}_{\star}$ (16)
$\displaystyle\dot{E}_{wind}=\alpha\dot{E}_{\star}$
$\displaystyle\dot{p}_{wind}=(\alpha\beta)^{1/2}\dot{p}_{\star}$

In the above equations $\beta$ accounts for ambient gas that is mixed into the
stellar ejecta (so that the total outflow rate in the wind fluid is a factor
$\beta$ larger than the rate at which stellar ejecta are created).555The
$\beta$ term is not to be confused with mass-loading factor
$\dot{M}_{out}/SFR$, where $\dot{M}_{out}$ is the total mass outflow rate in
all phases. For a standard IMF, $\beta=1$ implies that $\dot{M}_{\star}\sim
0.2$ SFR. The term $\alpha$ accounts for the effects of radiative losses that
drain away the energy carried by the wind fluid. Thus, $\beta\geq$ 1 and
$\alpha\leq$ 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 pre-existing gas clouds
that are being accelerated via momentum and/or kinetic energy transferred from
the wind fluid to the clouds 666There 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.. 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 $\alpha=\beta=1$, the 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 ($\tau_{1200}\sim 0.3$). This
means that only $\sim$ 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 $\alpha\sim\beta\sim 1$ (Strickland & Heckman
2009). In this case, $V_{wind}\sim$ 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
$\dot{E}_{\text{out}}$/$\dot{E}_{\star}$ are 1 to 20 %. We note that this does
not imply a small value for $\alpha$ (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 $\sim$ 30%. This suggests an efficient
transfer of momentum. Interestingly, there are several cases in which
$\dot{p}_{out}>$ $\dot{p}_{\star}$. This is not in conflict with the model of
an outflow driven by the wind fluid. That is, in the case where
$\alpha\times\beta>1$, the momentum flux in the wind fluid can exceed the
momentum flux from the starburst by a factor $(\alpha\times\beta)^{1/2}$ (see
Equation set (16)). This would correspond to a case in which radiative losses
in the wind fluid are negligible ($\alpha\sim$ 1) but the wind fluid is
significantly contaminated by ambient gas mixed into it ($\beta>>1$). This
would of course require an efficient transfer of this momentum from the wind
fluid to the clouds in the warm outflow.

###  6.2 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):

$\displaystyle\dot{p}_{crit}$ $\displaystyle=4\pi R_{0}N_{\text{H}}\mu
m_{p}v^{2}_{\text{cir}}$ (17)
$\displaystyle=10^{32.9}\text{dynes}\times\frac{R_{0}}{kpc}\times\frac{v_{\text{cir}}}{100kms^{-1}}\times\frac{N_{\text{H}}}{10^{20}cm^{-2}}$

The resulting $\dot{p}_{\text{crit}}$ values are listed in Table 5. In the
last panel of Figure 13, we compare the normalized outflow velocity, i.e.,
$V_{\text{out}}$/$V_{\text{cir}}$, versus the normalized outflow momentum
flux, i.e., $\dot{p}_{\text{out}}$/$\dot{p}_{\text{crit}}$. The two black
vertical lines split the figure into three different regimes for outflows: (1)
$\dot{p}_{\text{out}}$/$\dot{p}_{\text{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 $<$ $\dot{p}_{\text{out}}$/$\dot{p}_{\text{crit}}$ $<$ 10, where we
expect relatively weak outflows are driven in this regime. A few of our
observed outflows fall into this region. (3)
$\dot{p}_{\text{out}}$/$\dot{p}_{\text{crit}}$ $>$ 10, where relatively strong
outflows should be driven. We find most of our observed outflows fall into
this last region of $\dot{p}_{\text{out}}$/$\dot{p}_{\text{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_{\text{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_{\text{out}}$ = 100%, 50%, and 20%
of $V_{\text{max}}$, respectively. We find that most of our confirmed outflows
are consistent with the $V_{\text{out}}$ = 50% $V_{\text{max}}$ curve. This
can be explained as a “natural choice”, i.e., $V_{\text{max}}$ =
$V_{\text{out}}$ \+ 0.5 FWHM${}_{\text{out}}$ $\sim$ 2$V_{\text{out}}$ (see
Section 5.4.2). Unlike Heckman et al. (2015), we do not probe the regime of
low $\dot{p}_{\text{out}}$/$\dot{p}_{\text{crit}}$ very well, and so do not
see its correlation with $V_{\text{out}}$/$V_{\text{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.777We 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. In this
model, the gas we see in absorption is swept-up 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 $\dot{p}_{\text{out}}$/$\dot{p}_{\text{crit}}$ seen in Figure 13
above.

In this case, the predicted radius and expansion velocity are given by:

$\displaystyle R_{s}$
$\displaystyle=2.16~{}\dot{p}_{35}^{1/4}n_{0}^{-1/4}t_{7}^{1/2}kpc$ (18)
$\displaystyle=1.77~{}SFR_{10}^{1/4}n_{o}^{-1/4}t_{7}^{1/2}kpc$ $\displaystyle
V_{s}$ $\displaystyle=106~{}\dot{p}_{35}^{1/4}n_{0}^{-1/4}t_{7}^{-1/2}km/s$
$\displaystyle=87~{}SFR_{10}^{1/4}n_{0}^{-1/4}t_{7}^{-1/2}km/s$

Here $\dot{p}_{35}$ is the momentum flux in the wind in units of 1035 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 t7 is the time since the expansion began in
units of 107 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}\propto\dot{p}_{35}^{1/4}$ (i.e., $V_{s}\propto SFR^{1/4}$) is very
similar to the measured slope between $V_{\text{out}}$ and SFR ($\sim$ 0.24,
see Figure 9). In the second panel of Figure 9, we show the best-fit relation
with $n_{0}t_{7}^{2}$ treated as a free parameter (green dashed line). We find
a best-fit value of $n_{0}t_{7}^{2}$ = 2.8 $\times 10^{-2}$ cm-3 yr2. Assuming
$t_{7}\sim$ 1, this implies low-density gas (i.e., 2.8 $\times 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\times 10^{20}\dot{p}_{35}^{1/4}$
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 $\dot{p}_{*}$ 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\times 10^{-2}$ cm-3, we get:

$\displaystyle\dot{M}$
$\displaystyle=54~{}\dot{p}_{35}^{3/4}t_{7}^{1/2}M_{\odot}yr^{-1}$ (19)
$\displaystyle=31~{}SFR_{10}^{3/4}t_{7}^{1/2}M_{\odot}yr^{-1}$

These rates are about three times larger than we have estimated, and would
lead to a steeper slope in the dependence of $M_{\odot}$ 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${}_{\text{out}}$$<<$
$V_{\text{out}}$. This is inconsistent with the data (see Figures 2, 3, and
5).

###  6.3 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 sub-grid 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_{\text{out}}$ and $V_{\text{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_{\text{out}}$$\propto$
$V_{\text{cir}}^{0.6}$. 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 mass-loading factor
$\dot{M}_{out}/SFR$ should vary with $V_{\text{out}}$. Since these models
assume a linear relationship between $V_{\text{out}}$ and $V_{\text{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, $\dot{M}_{out}/SFR\propto
V_{out}^{-1}\propto V_{cir}^{-1}$ (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 $\dot{M}_{out}/SFR\propto
V_{out}^{-2}\propto V_{cir}^{-2}$. 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_{\text{out}}$ and $V_{\text{cir}}$ are both intermediate
between these two cases (slope $\sim$ –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_{\text{cir}}$ (providing better dynamic
range).

###  6.4 Comparisons to an Semi-Analytic Model of a Multi-phase Outflow

In this section, we compare our results to a new semi-analytic 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_{\text{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 $\beta$ for the wind fluid, which is
varied between 1 and 5, the initial mass-loading factor
($\dot{M}_{\text{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_{\text{out}}$ will be 260 km -1 based on the
SFR in the model and 360 km s-1 based on $V_{\text{cir}}$ in the model. These
velocities are most consistent with the most massive clouds in the models (106
M⊙), which reach $V_{\text{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 $\sim$ 0.6 based on SFR and $\sim$
0.3 based on $V_{\text{cir}}$. This is in the range adopted by the models for
the initial mass-loading factor. In the models, the mass-loading factor in the
warm outflow only drops significantly with radius for the case in which
$\beta=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 $\sim 30$ to 100% of $\dot{p}_{\star}$ for $\beta$
from $\sim$ 1 to 5, respectively. The $\beta=1$ model is therefore a better
fit to our data (Figure 13). Finally, the kinetic energy flux in the clouds
depends most strongly on $\beta$, being in the range 20 to 40%
$\dot{E}_{\star}$ for $\beta=5$ vs. only 4 to 10% for $\beta=1$. Our median
value of $\sim$5% (see Figure 13) again favors the $\beta=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 $\beta=5$ model, but only for the lower mass clouds in the $\beta=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 $\beta=1$ model
with massive clouds.

These models also make different predictions for the total hydrogen column
density ($N_{\mathrm{\scriptstyle H}}$) in clouds along the line of sight (D.
Fielding, private communication). For the $\beta=1$ models, $N_{H}$ increases
from about 2 $\times 10^{20}$ to 8 $\times 10^{20}$ cm-2 as the individual
cloud masses increase from $10^{2}$ to $10^{6}$ M⊙. For the $\beta=5$ model
the corresponding range is $\sim$ 1 to 2 $\times 10^{21}$ cm-2. For our sample
we find a median value of $N_{H}=5\times 10^{20}$ cm-2, which agrees best with
the $\beta=1$ model and massive clouds.

In summary, at least for parameters appropriate to an M82-level starburst, our
data are most consistent with the models with $\beta=1$ (a fast wind that has
not been contaminated) and with relatively massive clouds ($\sim 10^{5}$ to
106 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_{\text{cir}}$, SFR, and starburst size.

###  6.5 Comparisons to a Numerical Simulation of a Multi-phase Outflow

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_{\text{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\times 10^{5}$ K, 2) an intermediate phase with 5
$\times 10^{5}>T>2\times 10^{4}$ K, and 3) a cool phase with $T<2\times
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_{\text{out}}$ for the cool phase rise rapidly with
radius to a value of $\sim$ 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_{\text{out}}$ will be 260 km -1 based on their assumed SFR and 360 km s-1
based on their assumed $V_{\text{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\pi$ 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
mass-loading factor of 0.6 based on SFR and 0.3 based on $V_{\text{cir}}$.

Similarly, the momentum flux in the cool phase in the simulations is about 10
to 20% of $\dot{p}_{\star}$ (falling at large radii). This is a bit smaller
than the median value of $\sim$ 30% in the data (Figure 13). Finally, the
kinetic energy flux in the cool phase in the simulations is in the range of
$\sim$ 2 to 4% of $\dot{E}_{\star}$, 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 $\sim 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 ($\gtrsim$ a factor of 2) from Si ii to Si iv, which is not seen
in our galaxies (Figure 4). They also predict that $V_{\text{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_{\text{out}}$ correlates with the surface density of SFR with a
slope $\sim$ 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 $\sim$ 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.

Figure 14: Correlations between $V_{\text{out}}$ 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).

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.

###  6.6 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 ($\beta>>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).

Thompson et al. (2016) consider two cases. One is a case in which $\beta$ 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:

$\displaystyle V_{max}\sim 1250~{}km/s~{}[\alpha~{}\xi
SFR_{10}/((\Omega/4\pi)R_{*,0.3})]^{0.18}$ (20)

Here $\xi$ 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 $\beta$ is large enough for cooling to
occur at $R_{*}$.This yield a critical outflow velocity given by:

$\displaystyle V_{crit}\sim 720~{}km/s~{}[\alpha~{}\xi
SFR_{10}/((\Omega/4\pi)R_{*,0.3})]^{0.135}$ (21)

In Figure 15, we compare these two velocities to our observed values of
$V_{max}=V_{out}+0.5FWHM_{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 $\beta$ that the outflows would be still-born 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
velocities of the warm phase of $\sim 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.

Figure 16: Comparisons between the measured momentum flux of outflows
($\dot{p}_{\text{out}}$) to the predicted maximum wind momentum flux
($\dot{p}_{max}$) from Lochhaas et al. (2021). The three lines show Y = 10, 1,
0.1 X. There is only one galaxy that have $\dot{p}_{\text{out}}$ $>$
$\dot{p}_{max}$, which is similar to Figure 13. See Section 6.6 for
discussion.

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 ($\beta>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 $\dot{p}_{wind}=(\alpha\beta)^{1/2}\dot{p}_{*}$.
Lochhaas et al. (2021) calculated the maximum allowed value of the product of
$\alpha\beta$ under the conditions of increasingly large $\beta$ leading to
increasing strong radiative cooling ($\alpha<1$). They find the maximum wind
momentum flux to be:

$\displaystyle\dot{p}_{max}=7.2\times 10^{34}(R_{*,0.3})^{0.14}(\alpha
SFR_{10})^{0.86}dynes$ (22)

We compare our estimates of the outflow momentum flux with $\dot{p}_{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 $\dot{p}_{max}$ is generally similar to
$\dot{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 ($\alpha\beta\sim 1$).

##  7 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_{\text{out}}$) and the Full Width at Half Maximum (FWHM${}_{\text{out}}$)
of the blue-shifted absorption-lines, the total Si and H column densities
($N_{\mathrm{\scriptstyle Si}}$ and $N_{\mathrm{\scriptstyle H}}$), and the
outflow rates of Si ($\dot{M}_{\text{Si,out}}$), mass
($\dot{M}_{\text{out}}$), momentum ($\dot{p}_{\text{out}}$), and kinetic
energy ($\dot{E}_{\text{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\pi$ steradian (they cannot be
well-collimated).

  * •

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_{\text{out}}$ are quite consistent among all the
transitions we measure, there is significant scatter in the values for
FWHM${}_{\text{out}}$İn particular, we find a systematic trend for
FWHM${}_{\text{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_{\text{out}}$.

  * •

We found highly significant correlations of both $V_{\text{out}}$ and
FWHM${}_{\text{out}}$ with star-formation rate (SFR), galaxy stellar mass
(M⋆), and galaxy circular velocity ($V_{\text{cir}}$). The best-fit
relationship is $V_{\text{out}}$ $\propto$ $V_{\text{cir}}$0.6, and the ratio
of $V_{\text{out}}$/$V_{\text{cir}}$ being $\sim$ 6 and $\sim 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 $\sim$ 0.3). Most of the
observed extinction is produced by dust in the static ISM.

  * •

The average total hydrogen column density ($N_{\mathrm{\scriptstyle H}}$) is
$\sim 10^{20.7}$ cm-2, and that neutral hydrogen [N(H i)] only constitutes 0.1
to 1% of $N_{\mathrm{\scriptstyle 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_{\mathrm{\scriptstyle H}}$(v), we find the column densities peak
near $V_{\text{out}}$, while the broad wings of outflow profiles have
significantly smaller values of $N_{\mathrm{\scriptstyle H}}$.

  * •

We found a highly significant correlation between the mass outflow rate
($\dot{M}_{\text{out}}$) and the SFR, but with a shallow slope
($\dot{M}_{\text{out}}$$\propto SFR^{0.4}$). Hence the so-called mass-loading
factor of the outflow is $\dot{M}_{\text{out}}$/SFR $\propto SFR^{-0.6}$. We
also found that the mass-loading factor is a steep inverse function of both
$V_{\text{out}}$ and $V_{\text{cir}}$ (slope $\sim-1.6$), with mass-loading
factors of $\sim 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 ($\dot{p}_{\text{out}}$ and
$\dot{E}_{\text{out}}$) and the rates at which the starbursts supply momentum
and kinetic energy ($\dot{p}_{\star}$ and $\dot{E}_{\star}$). The median
values are $\dot{p}_{\text{out}}$ $\sim$ 30% $\dot{p}_{\star}$ and
$\dot{E}_{\text{out}}$ $\sim$ 5% $\dot{E}_{\star}$.

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-of-the-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 $\sim$ 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_{\text{out}}$ and $SFR$, it fails in other respects, and would produce
absorption-lines with $V_{\text{out}}$$>>$ FWHM${}_{\text{out}}$(while we find
$V_{\text{out}}$$\sim$ 0.5 FWHM${}_{\text{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 (105 to 106
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 $\sim$ 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,