Kinematics, Structure, and Mass Outflow Rates of Extreme Starburst Galactic Outflows

We present results on the properties of extreme gas outflows in massive (M * ∼ 1011 M ⊙), compact, starburst (star formation rate, SFR∼ 200 M ⊙ yr−1) galaxies at z = 0.4–0.7 with very high star formation surface densities (ΣSFR ∼ 2000 M ⊙ yr−1 kpc−2). Using optical Keck/HIRES spectroscopy of 14 HizEA starburst galaxies, we identify outflows with maximum velocities of 820–2860 km s−1. High-resolution spectroscopy allows us to measure precise column densities and covering fractions as a function of outflow velocity and characterize the kinematics and structure of the cool gas outflow phase (T ∼ 104 K). We find substantial variation in the absorption profiles, which likely reflects the complex morphology of inhomogeneously distributed, clumpy gas and the intricacy of the turbulent mixing layers between the cold and hot outflow phases. There is not a straightforward correlation between the bursts in the galaxies’ star formation histories and their wind absorption line profiles, as might naively be expected for starburst-driven winds. The lack of strong Mg ii absorption at the systemic velocity is likely an orientation effect, where the observations are down the axis of a blowout. We infer high mass outflow rates of ∼50–2200 M ⊙ yr−1, assuming a fiducial outflow size of 5 kpc, and mass loading factors of η ∼ 5 for most of the sample. While these values have high uncertainties, they suggest that starburst galaxies are capable of ejecting very large amounts of cool gas that will substantially impact their future evolution.


INTRODUCTION
In the last decade, the potential impact of galaxyscale outflows in galactic evolution has become widely recognized (e.g., Kormendy & Ho 2013;Somerville & Davé 2015;Veilleux et al. 2020). Outflows provide a mechanism that can regulate and possibly quench star formation activity in the galaxy by blowing away the gas that feeds star formation and supermassive black hole (SMBH) growth, enriching the large-scale galactic environment with metals. Theoretical studies and simulations show that this "feedback" offers a natural explanation for a variety of observations, e.g., the chemical enrichment of the circumgalactic and intergalactic medium (CGM, IGM), the self-regulation of the growth of the SMBH and of the galactic bulge, the relative dearth of both low and high mass galaxies in the stellar mass function, and the existence of the red sequence of massive, passive galaxies (e.g., Silk & Rees 1998;Thompson et al. 2005;Kereš et al. 2005;Oppenheimer et al. 2010;Murray et al. 2010;Tumlinson et al. 2017;Davies et al. 2019;Nelson et al. 2019). While galactic outflows appear to be crucial to rapidly shutting down star formation, the physical drivers of this ejective feedback are still unclear. In particular, the relative role of feedback from stars versus SMBHs in quenching star formation in massive galaxies remains widely debated.
Our team has been studying a sample of massive (M * ∼10 11 M ) galaxies at z = 0.4−0.8, originally selected from the Sloan Digital Sky Survey (SDSS; York et al. 2000) Data Release 8 (DR8; Aihara et al. 2011) to have typical signatures of young post-starburst galaxies. We refer to these galaxies as the HizEA 1 sample. Their spectra exhibit both strong stellar Balmer absorption from A-and B-stars and weak or absent nebular emission lines, implying minimal ongoing star formation. These galaxies are driving extremely fast ionized gas outflows, as revealed by highly blueshifted Mg II absorption in ∼90% (Davis et al., submitted) of their optical spectra, with velocities of ∼1000−2500 km s −1 , an order of magnitude larger than typical z ∼1 starforming galaxies (e.g., Weiner et al. 2009;Martin et al. 2012). These results initially led our team to conclude that these galaxies are post-starburst systems with powerful outflows that may have played a critical role in quenching their star formation (Tremonti et al. 2007).
However, several of these galaxies are detected in the Wide-field Infrared Survey Explorer (WISE; Wright et al. 2010), and their ultraviolet (UV) to near-IR spectral energy distributions (SEDs) indicate a high level of heavily obscured star formation (> 50 M yr −1 ; Diamond-Stanic et al. 2012). Furthermore, Hubble Space Telescope (HST) imaging of 29 of these galaxies reveals they are in the late stages of highly dissipational mergers with exceptionally compact central starforming regions (R e ∼few 100 pc; Diamond-Stanic et al. 2012Sell et al. 2014). Combining star formation rate (SFR) estimates from WISE rest-frame mid-IR luminosities with physical size estimates from HST imaging, we derive extraordinarily high SFR surface densities Σ SF R ∼10 3 M yr −1 kpc −2 (Diamond-Stanic et al. 2012), approaching the maximum allowed by stellar radiation pressure feedback models (Lehnert & Heckman 1996;Meurer et al. 1997;Murray et al. 2005; Thompson et al. 2005). These findings paint a different picture where our targets are starburst galaxies with dense, dusty star-forming cores (Perrotta et al. 2021), and a substantial portion of their gas and dust is blown away by powerful outflows. Millimeter obser-vations for two galaxies in our sample suggest that the reservoir of molecular gas is being consumed by the starburst with remarkable efficiency (Geach et al. 2013), and ejected in a spatially extended molecular outflow (Geach et al. 2014), leading to rapid gas depletion times. Notably, we find little evidence of ongoing AGN activity in these systems based on X-ray, IR, radio, and spectral line diagnostics (Sell et al. 2014;Perrotta et al. 2021). These results agree with models in which stellar feedback is the primary driver of the observed outflows (e.g. Hopkins et al. 2012Hopkins et al. , 2014. Our sample of compact starburst galaxies is characterized by the highest Σ SF R and the fastest outflows (>1000 km s −1 ) observed among star-forming galaxies at any redshift; therefore, they are an ideal laboratory to test the limits of stellar feedback and determine whether stellar processes alone can drive extreme outflows, without the need to invoke feedback from SMBH. They could represent a short but relatively common phase of massive galaxy evolution (Whalen et al. 2022).
Low and medium-resolution observations have shown that large-scale galactic outflows are a common feature of star-forming galaxies across a wide range of masses and redshifts (e.g., Martin 1998;Heckman et al. 2000;Weiner et al. 2009;Martin & Bouché 2009;Martin et al. 2012;Rubin et al. 2010Rubin et al. , 2014Heckman et al. 2015;Heckman & Borthakur 2016;McQuinn et al. 2019). Thanks to these studies we have a fair understanding of the typical outflow kinematics, and how their main properties scale with the fundamental properties of their galaxy hosts. While absorption-line spectroscopy at low and medium resolution is sufficient to measure velocities and the total amount of absorption within outflows, robust determinations of key physical properties such as the column density, covering fraction, and mass outflow rate require high spectral resolution. In this paper, we present new optical Keck/HIRES observations for 14 of the most well-studied starburst galaxies in our sample. These high-resolution spectra allow us to directly measure accurate column densities and covering fractions as a function of velocity for a suite of Mg and Fe absorption lines. We utilize this to probe the small scale structures of the extreme galactic outflows observed in our sample and investigate the potential impact of these outflows on the evolution of their host galaxies.
The paper is organized as follows: Section 2 describes the sample selection, observations, and data reduction; Section 3 illustrates our line profile fitting method and measurements of the absorption line column densities and kinematics; Section 4 presents our main results in comparison to theoretical models and other galaxy samples, and discusses the more comprehensive implications of our analysis. Our conclusions are reviewed in Section 5.

SAMPLE AND DATA REDUCTION
The parent sample for this analysis was selected from the Sloan Digital Sky Survey I (SDSS-I York et al. 2000) Data Release 8 (DR8; Aihara et al. 2011) and is described in Tremonti et al. (2007); Diamond-Stanic et al. (2012); Sell et al. (2014); Diamond-Stanic et al. (2021) and Davis et al., (submitted). Briefly, as mentioned above, the original goal was to identify a sample of intermediate redshift (z = 0.35−1) post-starburst galaxies to investigate their star formation quenching mechanisms. We focused on redshift > 0.35 so that the Mg IIλλ2796, 2803 doublet, an ISM line broadly used to probe galactic winds, was easily observable with optical spectrographs. Since the SDSS main galaxy sample (Strauss et al. 2002) is magnitude limited (r < 17.7) and does not contain many galaxies at z > 0.25, we concentrated on objects targeted for spectroscopy as quasars because they probe fainter magnitudes (i < 20.4) and higher redshifts. We selected all objects classified by the SDSS spectroscopic pipeline as galaxies with redshifts between z = 0.35 − 1 and either g < 20 or i < 19 mag (296 galaxies). To select galaxies with a recent (< 1 Gyr) burst of star formation, but little ongoing activity, we required them to show strong Balmer-absorption (suggestive of a burst 50 − 1000 Myr ago) and moderately weak nebular emission (suggestive of a low SFR in the past 10 Myr). Removing 3% of the objects by visual examination (blazars or odd quasars with wrong redshifts) yielded a sample of 121 galaxies. More details about the sample selection can be found in Tremonti et al., (in prep.) and Davis et al., (submitted).
We have undertaken extensive follow-up observations on a subsample of 50 of these galaxies. In choosing this subsample, we prioritized galaxies with the largest gband fluxes and the youngest stellar populations, but we included galaxies with older burst ages for comparison. The sample comprises galaxies with mean stellar ages in the range of 4 -400 Myr. Notably, we did not prioritize galaxies with Mg II absorption lines.
We obtained ground-based spectroscopy of the 50/121 galaxies with the MMT/Blue Channel, Magellan/MagE, Keck/LRIS, Keck/HIRES, and/or Keck/KCWI (Tremonti et al. 2007;Diamond-Stanic et al. 2012;Sell et al. 2014;Rupke et al. 2019), X-ray imaging with Chandra for 12/50 targets (Sell et al. 2014), radio continuum data with the NSF's Karl G. Jansky Very Large Array (JVLA/VLA) for 20/50 objects (Petter et al. 2020), millimeter data (ALMA) for 2/50 targets (Geach et al. 2014(Geach et al. , 2018, and optical imaging with HST for 29/50 galaxies ("HST sample" Diamond-Stanic et al. 2012;Sell et al. 2014;Diamond-Stanic et al. 2021). For the HST observations, we first targeted the twelve galaxies that were most likely to have AGN. We observed three galaxies with broad Mg II emission (expected to be type I AGN) and nine galaxies with strong [O III] emission, which could be indicative of an obscured (type II) AGN. We subsequently observed an additional seventeen galaxies that have the youngest estimated post-burst ages (t burst < 300 Myr). We also obtained multi-band HST imaging to study the physical conditions at the centers of the 12/29 galaxies with the largest SFR surface densities estimated by Diamond-Stanic et al. (2012), (30 M yr −1 kpc −2 < Σ SF R < 2000 M yr −1 kpc −2 ), and investigated the young compact starburst component that makes them so extreme (Diamond-Stanic et al. 2021).
In this paper, we focus on 14 galaxies from the HST sample (6 from the 12/29 most AGN-like galaxies, and 8 from the 17/29 with the youngest post-burst ages). The galaxies and their basic properties are listed in Table 1.

Data Acquisition and Reduction
We obtained high-resolution spectra of 14 starburst galaxies spanning emission redshifts 0.4 < z < 0.8, using the Keck High Resolution Echelle Spectrometer (HIRES; Vogt et al. 1994) on the Keck I telescope over the nights of 12 − 13 August 2009, 10 − 11 December 2009, and 13 − 14 May 2012. HIRES was configured using the 1.148 × 7 arcsec C5 decker with no filters and 2x2 binning. The 1.148 arcsec slit width disperses the light with R = λ/∆λ = 37,000 ( 8 km s −1 ) projected to 3 pixels per resolution element, ∆λ. Individual exposures were 3600 seconds, with total integration times of 2 − 3 hours per object.
This yielded an observed wavelength coverage from ∼3150 to 5800Å in >30 orders, which for our z = 0.4 − 0.8 galaxies corresponds to a rest-frame coverage from λ rest ∼1800−3250 to 2300 − 4100Å depending on redshift. We concentrate our analysis on the Mg IIλλ2796, 2803 doublet, the Mg Iλ2852 transition, and the Fe IIλ2344, 2374, 2582, 2587, and 2600 multiplet. The data were reduced using the XIDL HIRES Redux code 2 . We used this reduction pipeline to perform flat-field corrections, determine wavelength solutions using thorium-argon (ThAr) arcs, measure slit profiles for each order, perform sky subtraction, and remove cosmic rays. The output of the spectral extraction is an orderby-order spectrum of counts versus vacuum wavelength. For a more comprehensive description of this procedure, we refer the reader to O'Meara et al. (2015). As automated flux calibration of HIRES spectra is unreliable (e.g., Suzuki et al. 2005), each echelle order was examined by eye to perform continuum fitting. A Legendre polynomial was fit to each order and adjusted using anchor points in absorption-free regions with a custom IDL code.

Galaxy properties
We report in Table 1 relevant derived galaxy properties for our sample. Accurate systemic redshifts are nec-2 Available at https://www.ucolick.org/ ∼ xavier/HIRedux/index. html essary in order to determine outflow velocities. For this purpose, we use stellar continuum fits as derived in Davis et al. (submitted). We have collected low/medium resolution (R = 600 − 4000) optical spectra of our parent sample (50 galaxies) using three instruments on 6 − 10-m class telescopes (MMT/Blue Channel, Magellan/MaGE, and Keck/LRIS; see Davis et al. (submitted) and Tremonti et al. in prep. for details on these data sets). We required our spectra to extend to at least 4000Å in the restframe to cover strong stellar or nebular features to aid in measuring the galaxy systemic redshift (e.g. [O II], low order Balmer absorption lines). We fit the spectra with a combination of SSP models and a Salim et al. (2018) reddening law derived from local galaxies that are analogs of massive, high redshift galaxies. We used the Flexible Stellar Population Synthesis code (Conroy et al. 2009;Conroy & Gunn 2010) to generate SSPs with Padova 2008 isochrones (Marigo et al. 2008), a Salpeter (1955 initial mass function (IMF), and the theoretical stellar library C3K (Conroy et al., in prep.) with a resolution of R ∼10,000. We utilize solar metallicity SSP templates with 43 ages spanning 1 Myr − 8.9 Gyr and perform the fit with the Penalized Pixel- Fitting (pPXF) software (Cappellari & Emsellem 2004;Cappellari 2017). We mask the region around Mg II and the forbidden emission lines during the fit and implement two separate templates for broad and narrow Balmer emission lines, assuming Case B recombination line ratios. Both the lines and stellar continuum are attenuated by the same amount of dust in the pPXF fit; we find that the results are insensitive to this assumption. By fitting Balmer emission and absorption lines simultaneously we can take into account the potential infill of the absorption line cores. One of the outputs of our pPXF fit is the galaxy systemic redshift (z sys ), which we list in Table 1. The fit produces a stellar continuum model without a nebular emission component (see Figure 1). Most sources, in addition to having strong Balmer absorption, have very blue continua indicating a recent starburst event (∼1−10 Myr) that is not highly dust-obscured. For galaxies with young stellar populations, like those in our sample, the stellar contribution to the Mg absorption lines is minimal. However, we utilize our best fit pPXF continuum model to properly remove the stellar absorption features from each spectrum in our sample.
To estimate stellar masses (M * ) and star formation rates (SFR) we fit the combined broad-band UV -mid-IR photometry and optical spectra using the Bayesian SED code Prospector (Leja et al. 2019;Johnson et al. 2021), as described in Davis et al. (submitted). Briefly, we incorporate the 3500 -4200Å spectral region in the SED fit as it covers many age-sensitive features (e.g., D4000, Hδ). We generate simple stellar population (SSP) models employing the Flexible Stellar Populations Synthesis code (FSPS; Conroy et al. 2009), adopting a Kroupa IMF (Kroupa 2001) and the MIST isochrones (Choi et al. 2016) and the C3K stellar theoretical libraries (Conroy et al., in prep.). The stellar models are very similar to the ones described above over the wavelength range of interest for this work. We determined the best-fit parameters and their errors from the 16th, 50th, and 84th percentiles of the marginalized probability distribution function. The combined photometry and spectra are well fit by these models (see Davis et al., submitted for examples of the SED fitting). However, the dust emission properties of our sample are poorly constrained due to the low signal-to-noise ratio (SNR) of the WISE W3 and W4 photometry and the limited infrared coverage of the SED. This yields fairly tight constraints on the M * (±0.15 dex) and slightly larger errors on the SFR (±0.2 dex). M * represents the present-day stellar mass of the galaxy (after accounting for stellar evolution) and not the integral of the star formation history (i.e. total mass formed). We list in Table 1 the SFRs calculated from the resulting star formation history (SFHs; see Fig. 9), averaged over the last 100 Myr. This is the typical timescale for which UV and IR star formation indicators are sensitive (Kennicutt & Evans 2012). As here we are interested in the most recent SFH, we compute the light-weighted age of the stellar populations younger than 1 Gyr, with the light contribution estimated at 5500Å. These <1 Gyr lightweighted ages more closely replicate the timescale of the peak SFR than the mass-weighted ages. We use lightweighted ages rather than the age since the burst as the light-weighted ages are more robust to changes in our modeling approach. However, there may be systematic errors related to the stellar population models we adopt. In practice uncertainties in the treatment of Wolf-Rayet stars and high mass binary evolution can have a large impact on the UV spectra of galaxies with young stellar populations (e.g. Eldridge & Stanway 2016). The detailed analysis required to produce a quantitative estimate of the systematic errors on the light-weighted ages is beyond the scope of this work. The light-weighted ages for the galaxies in our sample are reported in Table 1.
The effective radii (r e ) measurements for the galaxies in our sample are discussed in Diamond-Stanic et al. (2012. Briefly, for 11 galaxies we use the GAL-FITM software Vika et al. 2013) to perform Sérsic fits of joint multi-band HST imaging (Diamond-Stanic et al. 2021) in the UVIS/F475W and UVIS/814W filters. To prevent uncertainties due to tidal features, we fit the central region of the galaxy and extrapolate the fit to larger radii to compute r e . The shorter-wavelength filter (λ rest (F475W) ≈ 3000Å) traces the young, unobscured stars, while the longerwavelength filter (λ rest (F814W) ≈ 5200Å) is more sensitive to the underlying stellar mass. Characteristic errors on the effective radius are of the order of 20%. For the remaining three galaxies (J1125, J1232, and J1450) we quantify the morphology utilizing optical HST UVIS/F814W image (Diamond-Stanic et al. 2012). We model the two-dimensional surface brightness profile with a single Sersic component (defined by Sersic index n = 4 and r e ) using GALFIT (Peng et al. 2002(Peng et al. , 2010. We adopt an empirical model point-spread function (PSF) produced using moderately bright stars in our science images.

SPECTRAL ANALYSIS
In this section, we briefly present our data, discuss the method and assumptions adopted for the line profile fitting and describe the techniques used to estimate errors. The line fitting results for each galaxy are listed in Tables 2 and 3.

Absorption Lines as Tracers of Galactic Outflows
Galactic winds are typically recognized through their kinematic signatures. Winds seen in absorption are identified as troughs detected in the foreground of the galaxy stellar continuum, blueshifted with respect to the galaxy systemic velocity (e.g., Martin & Bouché 2009;Martin et al. 2012;Kornei et al. 2012;Rubin et al. 2014;Prusinski et al. 2021). In this context the line velocity shift relative to the systemic redshift (z sys ) at which 98% (v 98 ) or 50% (v avg ) of the equivalent width (EW) accumulates moving from red (positive velocities) to blue (negative velocities) across the line profile can be used to parameterize the kinematics of absorption lines.
The galaxies in our sample commonly display extreme kinematics in the Mg II absorption lines. Figure 2 shows the continuum-normalized Mg II spectral region for each galaxy and illustrates the range of spectral absorption properties in our data. The Mg IIλλ2796, 2803 absorption lines are produced by the blending of multiple velocity components. Additionally, the Mg IIλ2796 line has blueshifted components that blend with the Mg IIλ2803 components in almost all of the galaxies (with the exception of J1125 and J2116). The velocity blueshifts that cause this blending range from several hundred up to over 2,000 km s −1 in the Mg II absorbing gas. Given the complex kinematics, we cannot directly integrate the spectrum to measure the EW of the single absorption components. However, we can calculate the total Mg II EW by integrating the normalized spectrum over the velocity range for which we detect Mg II absorption. The total rest frame EWs for the galaxies in our sample are reported at the bottom of each panel in Figure 2. The Mg II absorption lines are very strong, with EWs from 1.89 to 12.90Å, with an average value of 7.49Å. (Below in Section 4.4 we compare our galaxies to other samples in the literature.) The v 98 values in the Mg II absorption lines span a range from −620 to −2700 km s −1 , with an average value of −1630. Such large line blueshifts are unambiguous signs of outflowing gas. We note that the potential emission filling effects on the v 98 estimates are (v = 0 km s −1 ) and Mg IIλ2803 (v = +770 km s −1 ) position in this velocity space, at the galaxy systemic redshift. We report the total rest frame Mg II equivalent width at the bottom of each panel. The spectra here and below have been lightly smoothed by two resolution elements (∼16 km s −1 ).
negligible (see Section 4.2). In the following section, we construct a model that describes our data in order to quantify the kinematics and strength of the outflows.

Line Fitting
We refer the reader to Rupke et al. (2005a) for a complete discussion of the analysis techniques used in absorption line fitting, depending on the type of line profiles studied. Here we summarize the relevant aspects for our sources which have partially covered, blended absorption lines.
The profile of a single absorption trough is parameterized by the distribution of the optical depth (τ ) along the line of sight, combined with how the absorbing gas covers the background source. If it completely covers the background source, such that the covering fraction (C f ) is unity, then the τ distribution can be calculated accurately as a function of velocity. This is valid also for a blended doublet, triplet, or higher-order multiplet lines which can be fit by solving a set of linear equations, using a "regularization" method (Arav et al. 1999). However, when the background light source subtends a wide angle on the sky compared to the absorbing clouds, C f is generally less than unity. This latter case may apply to our sample, where the central starburst illuminates gas clouds in the galaxy halo.
The intensity of an absorption line, where C f and τ depend on velocity and the continuum level has been normalized to unity, can be described as I(λ) = 1 − C(λ) + C(λ)e −τ (λ) . There is a degeneracy when solving for τ and C f for a single line. This degeneracy can usually be resolved by simultaneously fitting two or more transitions of the same species with known oscillator strengths (f 0 ), as the relative depths of the lines are independent of C f . However, in the case of doublet or higher-order multiplet lines blended together, we cannot solve for τ and C f directly, as in the region of overlap the solution is not unique. Nor can we directly integrate the spectrum to estimate the EW of the single absorption troughs. We must thus fit analytic functions to the line profiles. We, therefore, use intensity profiles which are direct functions of physical parameters (i.e., velocity, optical depth, and covering fraction). We assume τ can be formulated as a Gaussian, τ (λ) = τ 0 e (λ−λ0) 2 /(λ0b/c) 2 , where λ 0 and τ 0 are the central wavelength and central optical depth of the line, and b is the velocity width of the line (or Doppler param- ln 2]). The clear advantage of this approach is that the derived profile shapes are readily interpreted in terms of these physical parameters. This method accurately handles the intensity profile at both low and high optical depth in the case of constant C f which is an assumption we make for simplicity.
To decompose the blended absorption profiles, we need to make assumptions about the geometry of the absorbing gas. Here we refer to "components" to describe distinct velocity components of the same transition.
We assume the case of completely overlapping atoms when combining the intensities of two doublet or multiplet lines within a given velocity component. In this context, the atoms at all velocities are placed at the same position in the plane of the sky relative to the background continuum source. This is a simplification, as the broad profiles we observe could be due to the large-scale motions of individual clouds that are not all coincident. The covering fraction in this case is independent of velocity. The expression for the combined intensity of a doublet or multiplet (each with optical depth τ i (λ) and covering fraction C f,i (λ)) is then given by (1) We assume the case of partially overlapping atoms when combining two different velocity components. The motivation for this is that different components could have different C f if we assume they are spatially distinct. In this scenario, at a given wavelength there is an overlap between the atoms producing the components, and the covering fraction describes the fractional coverage of both the continuum source and the atoms producing the other component. The definition of the intensity is We fit our data utilizing a suite of IDL routines from the IFSFIT library (Rupke 2014). The code combines the two cases illustrated above when simultaneously fit-ting doublet lines that have multiple blended velocity components. We describe the normalized line profile intensity for each component as a function of four parameters: λ 0 , b, C f , and column density (N ). N can be formulated as Different transitions of the same ionic species are required to have the same component structure and are fit simultaneously to produce a single solution. However, we do not impose the same kinematics on different ionization states of the same element or on other elements. As we impose no constraints on the decomposition between different ionization states or different ionic species, any qualitative or quantitative correspondences between the fits occur naturally.
A comparison of the doublet (multiplet) line shapes in our spectra often reveals nearly identical intensities. In these cases, the relative intensities of the doublet (multiplet) troughs constrain that transition to be optically thick, setting a lower limit on the optical depth at the line center. We also enforce an upper limit on the optical depth in Mg II of τ 2803 = 10 and in Fe II of τ 2586 = 10. The data are not sensitive to changes in optical depth above these values unless the optical depth is 10. Our best fit is degenerate with fits to larger numbers of components. For example, an absorption line with a single component could be fit with a given b and τ 0 , or it could be fit by adding multiple velocity components with narrower b and larger τ 0 . Following previous studies (Rupke et al. 2005a;Martin & Bouché 2009), here we adopt the minimum number of velocity components required to describe the doublet (multiplet) absorption troughs.

Mg II
To quantify the kinematics and absorption strength of Mg II lines in our spectra, we fit the line profile shape assuming that the absorption of the continuum emission is due to foreground Mg II ions. However, continuum photons can be absorbed by gas both in front of and surrounding the galaxy and successively re-emitted in any direction, such that the excited ions decay straight back to the ground state. This scattering mechanism can produce P Cygni-like line profiles for Mg II Prochaska et al. 2011) which has often been detected in star-forming galaxies at z ∼0.3 and ∼1 (Martin & Bouché 2009;Weiner et al. 2009;Rubin et al. 2010Rubin et al. , 2014Rupke et al. 2019;Burchett et al. 2021). We observe Mg II emission in 9/14 galaxies in our sample. For 7 of these galaxies, we have KCWI data that confirm the presence of Mg II emission (see Section 4.2). Another common characteristic in our sample is the absence of strong Mg II absorption at the galaxy systemic velocity, which traces absorption in the interstellar medium (ISM) of the galaxy. We note that 11/14 galaxies in our sample have less than four percent of the Mg II EW within 200 km s −1 of z sys . We come back to this point and discuss the potential effects of Mg II emission filling below in Section 4.2.
We quantify the Mg II kinematics for each galaxy in our sample by fitting the absorption line profiles as described in Section 3.2. At a given wavelength the atoms producing the Mg IIλ2796 line are separated by 7Å (770 km s −1 ) from the ones producing the Mg IIλ2803 line. However, at a given velocity relative to z sys , we consider them to have (1) relative τ defined by atomic physics (τ 2976 = 2 × τ 2803 ; Morton 2003), and (2) equal b and C f . The best fit parameters are reported in Table 2. Figure 3 shows an example of the best fit to the Mg II absorption in the galaxy J1232+0723 (hereafter J1232). The lower horizontal axis shows the velocity of the Mg IIλ2796 doublet component relative to the galaxy systemic redshift (z sys ), and the vertical black dotted lines mark the Mg IIλ2796 and Mg IIλ2803 location in the velocity space at z sys , respectively. The Mg II doublet absorption profile is well described by six Gaussian components (shown with light blue solid lines). The brackets in Figure 3 mark Mg II doublet pairs numbered in decreasing order of blueshift velocity from z sys . The three most blueshifted components (1, 2, and 3) are characterized by C f < 1. Comparison of the Mg IIλ2796 and Mg IIλ2803 trough shapes illustrates their nearly identical intensity, but the lines are not black. The remaining three components have C f = 1. The spectrum also clearly includes redshifted (∼410 km s −1 ) Mg IIλ2803 emission observed up to +900 km s −1 relative to z sys . (The corresponding Mg IIλ2796 line is not obvious because of Mg IIλ2803 absorption at the same wavelengths.) The inclusion of an emission component in the model significantly improves the fit to the Mg II absorption trough. The red dashed line in Figure 3 represents the total Mg II best fit, including the Mg II emission (dark blue solid line). Figure 4 shows an ionic species stack from the spectrum of J1232, centered on the absorption lines included here. Each panel shows the line decomposition (light blue solid lines) and the total best fit (dashed orange lines) to a different ion absorption trough. The light grey part of the spectrum shown in some panels marks a region where multiple transitions overlap in wavelength. The top two left panels show our best fits to Mg IIλ2796 and Mg IIλ2803, respectively. It is worth noting that where the Mg IIλ2796 and Mg IIλ2803 velocity components lie on top of each other, then the dashed orange line does not describe the observed spectrum (in black) because it represents the total fit of a single transition only. For example, the high-velocity Mg IIλ2803 absorption makes the Mg IIλ2796 trough appear to be deeper than the Mg IIλ2803 trough between +190 and +330 km s −1 . When the best fit profile of both doublet components are combined they fit the data well, as shown in Figure 3 (dashed red line).

Mg I
We also fit Mg I absorption in our sample. Observations of a single Mg I transition do not directly constrain the Mg Iλ2853 optical depth. To quantify the kinematics and absorption strength of Mg Iλ2853 line profiles in our spectra, therefore, we follow two distinct approaches. First, we perform a fit to the Mg I absorp-tion troughs independent from the Mg II fit, as described in Section 3.2. We set C f = 1 and let the other parameters (i.e. b, λ 0 , and τ 0 ) vary. This approach assumes that an absorption line that is not black is produced by optically thin gas. If the gas is instead optically thick and the shape of the absorption trough is determined by the C f this method provides a lower limit to its column density. The Mg I best fit parameters are listed in Table 3 (columns 5, 6, and 7).
Our second approach assumes the optical depth of each Mg I velocity component to be linked to that of the Mg II component. We set b and λ 0 for each Mg I velocity component to be the same as the corresponding Mg II component and estimate what C f (Mg I) best fits the data. The optical depth in Mg I is (4) where χ(Mg I)/χ(Mg II) is the relative ionization correction. For starburst spectral energy distributions (SED), χ(Mg II) ≥ 0.7 ). We set the neutral fraction to be 30% (i.e. χ(Mg I) = 0.3). This approach is motivated by the fact that the Mg I and Mg II absorption troughs have similar kinematic structures, however, Mg I is seen to be shallower than Mg II. When τ 0 (Mg II) is large and C f (Mg II) < 1, the neutral Mg fraction must be less than a few percent to produce a shallow Mg Iλ2853 trough with C f = 1 (i.e. optically thin). In this case, it is more likely that Mg Iλ2853 is optically thick, and the C f determines the shape of the Mg I absorption trough. We report the results of this approach in Table 3 (columns 2, 3, and 4). We refer to the first approach to the Mg I fit as "independent", and to the second as "constrained". Figure 4 shows an example of our best fit to the Mg I lines in the galaxy J1232. The third from the top left panel shows the best fit to Mg Iλ2853. Following the independent approach, we identify three components falling within −380 and +220 km s −1 of z sys . They have good kinematic correspondence to three Mg II components, but the Mg I profiles are shallower. In Figure 4 (and for the rest of the sample Figure 5 − 24) we show the Mg Iλ2853 best fit as derived following the constrained approach. Even in this case, we find that three Mg I velocity components within −380 and +220 km s −1 of z sys provide a good fit to the data. The three Mg I components are described by C f values in the range of 16 − 33% of C f (Mg II).

Fe II
The spectral coverage of our data includes a series of strong Fe II resonance lines which have some useful advantages over analysis of the Mg II doublet. First, the large number of Fe II transitions (i.e. Fe IIλ2344Å, Fe IIλ2374Å, Fe IIλ2382Å, Fe IIλ2586Å, and Fe IIλ2600Å) span a wide range of oscillator strengths, which makes it possible to place robust bounds on the column density of singly-ionized iron, thereby better constraining the total gas column density. For example, the Fe IIλ2374 oscillator strength is a tenth that of the strongest transition, Fe IIλ2382Å (which is nearly equal to that of Mg IIλ2803Å). Furthermore, the absorption of Fe II in several of these transitions is followed by fluorescence (rather than resonance) emission, providing a clearer view of the intrinsic absorption profile Erb et al. 2012;Martin et al. 2012). The Fe IIλ2374Å transition is particularly useful as the resonance absorption trough is not filled in by the emission of fluorescent Fe II * λ2396Å photons.
We quantify the Fe II kinematics for each galaxy in our sample by fitting the absorption line profiles as described in Section 3.2. We model the observed transitions simultaneously to constrain C f and τ 0 . At a given velocity relative to z sys , we consider the transitions to have (1) a relative τ defined by atomic physics (Morton 2003), and (2) equal b and C f across the transitions. We do not detect any Fe II emission. The Fe II best fit parameters are reported in Table 2. In some spectra we additionally detect and model the Mn II λλλ 2576, 2594, and 2606 triplet which can blend with the Fe IIλλ 2586, 2600 transitions. The Mn II best fit parameters are also reported in Table 4.  As Fe II provides a more reliable estimate of the central optical depth, for the velocity components that are detected in both Fe II and Mg II we compare τ 0 (Mg II) derived from the Mg II fit (as described in Section 3.2.1) with the optical depth inferred from the Fe II fit. We assume a solar abundance ratio, which is supported by an ensemble of line ratio diagnostic diagrams used to estimate the metallicities in our sample (Perrotta et al. 2021). For a solar abundance ratio (log Mg/H = −4.45; Asplund et al. 2021), and conservatively adopting a comparable ionization correction (e.g. χ(Mg II) ≈ χ(Fe II)), the optical depth in the weaker magnesium doublet line is 10 −4.45 10 logMg/H × 10 logFe/H 10 −4.54 While the relative optical depth in Fe IIλ2382 and Mg IIλ2803 can be similar to each other, the Mg II optical depth can be much larger as Fe is more depleted onto dust grains than Mg (Savage & Sembach 1996;Jenkins 2009). Since there are not many studies on dust depletion in galactic winds, we adopt dust depletion factors (d(X)) consistent with those measured in the Galactic ISM (0.5 dex for Mg and 1.0 dex for Fe; Jenkins 2009). We note that Jones et al. (2018) studied a sample of nine gravitationally lensed z ≈ 2 − 3 star-forming galaxies, and inferred that the outflowing medium is characterized by moderate dust depletion d(Fe) = −0.9 dex, in line with our adopted d(Fe). Assuming a similar dust depletion correction for Mg and Fe, we would infer a τ 0 (Mg II 2803 ) that is systematically lower by a factor of three; this would decrease the inferred Mg II column densities (presented below) by roughly a factor of three. Adopting the dust depletion typical of the Galactic halo (0.59 dex for Mg and 0.69 dex for Fe; Savage & Sembach 1996) would lead to a similar result, with inferred Mg II column densities lower by a factor of 2.3. Hereafter we refer to the Mg IIλ2803 central optical depth derived directly from the Mg II fit as "measured" (τ 0 (Mg II) meas ), and the optical depth deduced using Equation 5 as "inferred" (τ 0 (Mg II) inf ). Figure 4 shows an example of our best fit to the Fe II lines in the galaxy J1232. The bottom left and the right column panels show the best fit to the Fe II transitions observed in this study. The shape of the Fe II absorption trough is well described by three kinematic components within −380 and +220 km s −1 of z sys . Although we do not force agreement in the kinematics of Fe II and Mg II, the different velocity components agree remarkably well (with comparable v and b to within the errors) over the velocity range where we detect Fe II. We can therefore directly compare their C f . The three Fe II components all have lower C f with values in the range of 46 − 71% of C f (Mg II). The τ 0 (Mg II) meas is not consistent with τ 0 (Mg II) inf for all the velocity components. From this comparison, we infer that the actual Mg II optical depth is ∼4 times larger than the minimum required to fit the doublet. The Mg II velocity components detected at −67 and −247 km s −1 are saturated and the best fit τ 0 values represent lower limits. For the component detected at 128 km s −1 , we ascribe the difference as resulting from uncertainties in the emission filling affecting this component.

Individual Systems
In the sections below we present fits to the Mg II, Mg I, and Fe II absorption lines detected in other three galaxies in our sample, to illustrate the range of line properties observed in our data. Fits for the remaining galaxies are presented in Appendix A. We generally find that the kinematics of Fe II match those of Mg II well, for the components with high EW and SNR.
3.3.1. J0905+5759 Figure 5 shows the best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J0905+5759 (hereafter J0905). We present the absorption line fit for this source as an example that exhibits relatively simple kinematics, with few components, and with Mg IIλλ2796, 2803 lines that do not blend together. There is strong agreement in the component structure of all the absorption lines. The shape of the Mg II, Mg I, and Fe II troughs are well described by two velocity components spanning −2900 to −2000 km s −1 of z sys . The spectrum also displays slightly redshifted (+56 km s −1 ) Mg IIλλ2796, 2803emission spanning −300 to +400 km s −1 . However, the resonance emission does not fill in the high-velocity Mg II absorption lines. The comparable intensity at all velocities of the Mg IIλ2796 and Mg IIλ2803 troughs indicates the Mg II is optically thick. As the lines are not black (though the lower velocity component is close to black) the covering fraction (C f < 1) determines the shape of the absorption troughs.
The independent Mg I fit identifies two components with similar kinematics to the constrained fit within the errors. As seen in J1232 above, here the Mg II troughs are deeper than the Mg I troughs, which indicates that only a fraction of the cloud volume contains much neutral Mg. It seems more likely that the Mg I absorption profile shape is determined by C f rather than optical depth, since the latter (i.e. C f (Mg I) = 1) would require the Mg neutral fraction to be less than a few percent (see Formula 4). The covering fraction for Mg I is 67 and 40% that of the corresponding Mg II velocity components. Before performing the Fe II fit, we identify and model the Mn II λλλ 2576, 2594, and 2606 triplet (yellow solid lines). In this case, they do not blend with the Fe IIλλ 2586, 2600 transitions. The Fe II transitions show roughly equal intensities, suggesting they are also tracing optically thick gas. The two Fe II velocity components are not black, and we find they are well described by C f (Fe II)< 1 (90% and 40% of C f (Mg II), respectively). Based on a comparison of τ 0 (Mg II) meas and τ 0 (Mg II) inf , we conclude that the Mg II fit provides only a lower limit on N(Mg II) and that the actual central optical depth for the two Mg IIλ2803 components is closer to 79 and 67 (rather than 6 and 4 as measured from the Mg II fit). Figure 6 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J1450+4621 (hereafter J1450). This spectrum exhibits a strong and complex Mg II absorption trough within −2000 and +350 km s −1 of z sys . The Mg IIλ2796 and Mg IIλ2803 line profiles are well described by ten velocity components. Many components (seven in Mg IIλ2796, and five in Mg IIλ2803) blend together. We do not detect Mg II emission in this galaxy. We find that two Mg II components (at v = +198 and v = −50 km s −1 ) are optically thick, and their shape is therefore determined by the covering fraction (C f = 0.15 and C f = 0.95). The remaining eight blueshifted components trace optically thin gas with C f = 1. We detect a strong Fe II absorption line at v = −45 km s −1 of z sys , that has good kinematic correspondence to the deepest Mg II component. We model this using only one component and find b(Fe II) to be 26% narrower than b(Mg II). The absorption trough is black for all Fe II transitions other than Fe IIλ2374. This provides a robust constraint on N(Fe II). Comparing τ 0 (Mg II) meas and τ 0 (Mg II) inf for the only component detected in both Mg II and Fe II, we infer that the Mg II fit provides a lower limit on N(Mg II), and the real τ 0 (Mg IIλ2803) is closer to 74 (than 8, as measured from the Mg II fit). We identify one Mg I absorption trough close to z sys (at v = −38 km s −1 ) that is visibly less deep than the corresponding Mg II and Fe II components. The independent and constrained fits agree well in terms of the absorption line kinematics. However, as τ 0 (Mg II) is at least ∼8 (and more likely 19), we conclude that the constrained fit provides only an upper limit to N(Mg I) and that the Mg I absorption is tracing optically thick gas (with C f (Mg I) = 32% C f (Mg II)). Mg I shows a second weak absorption trough at v ∼-1700 km s −1 that has a good kinematic correspondence to the two most blueshifted Mg II components. We model this using two velocity components and find that the independent fit results in a broader profile than that of the constrained fit. However, the kinematics agree within the errors given the large uncertainties in the independent fit, due to the poor SNR. As the Mg II in these two components traces optically thin gas, we conservatively favor the Mg I independent fit over the constrained fit. Figure 7 shows the best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J1341-0321 (hereafter J1341). This spectrum includes Mg II redshifted emission (+355 km s −1 ) within −200 and +1000 km s −1 of z sys . We detect strong absorption from Mg II falling within −2300 and +250 km s −1 of z sys . We model the Mg IIλ2796 and Mg IIλ2803 troughs using nine velocity components that are characterized by complex kinematics. Most of the components (eight Mg IIλ2796, and six Mg IIλ2803) blend together. The Mg II best fit includes a combination of optically thin (in three components) and optically thick (in five components) gas, and all components have C f < 1. Of note, we identify three extremely narrow redshifted Mg II components. These lines lie on top of the Mg II emission and are not resolved; their best fit b parameter is similar to one resolution element (∼8 km s −1 ). We detect Mg I and Fe II absorption within −1500 and −200 km s −1 of z sys , which we model using four and three velocity components, respectively. There is strong agreement in the Mg I kinematics of the independent and constrained fit solutions (with comparable v and b within the errors). The inde-

J1341-0321
There is a remarkable correspondence between Mg II and Fe II. However, the most blueshifted Fe II component is ∼57% broader than the corresponding Mg II. We ascribe the difference to the lower SNR in the Fe II spectral region compared to Mg II, which results in fitting Fe II using one component instead of two as for Mg II. We find Fe II to have lower C f than Mg II (24 − 77 % of C f (Mg II)). Based on a comparison of τ 0 (Mg II) meas and τ 0 (Mg II) inf , we conclude that the Mg II fit provides only a lower limit on N(Mg II) and that the actual central optical depth for the three Mg IIλ2803 components is closer to 17, 21, and 12 (rather than 4, 1, and 5 as measured from the Mg II fit). We detect outflows in 14 compact, massive starburst galaxies in absorption from Mg I, Mg II, and Fe II and show remarkably similar profiles of the absorption troughs in all these transitions. In most of our sample, the velocity dependence of the gas covering fraction (C f ) across components defines the absorption trough profile. Similarities in the profiles suggest these species reside in the same low-ionization gas structures. However, Mg II has on average a higher C f than Fe II at a given velocity, and a higher C f than neutral Mg, implying that the absorbing clouds or filaments are not homogeneous.
We now discuss our results, including the variation of the absorption line profiles and the possible connection with the star formation history of each galaxy (Section 4.1). In Section 4.2 we examine the lack of substantial absorption at the systemic redshift and the potential effect of emission line filling. We then use our results to gain insights into the role of different physical mechanisms in the extreme galactic outflows observed in our sample (Section 4.3). In Section 4.4 we investigate trends (or lack thereof) between the outflow absorption strength and galaxy properties. We conclude by presenting estimates of the mass outflow rates in our sam-ple and discussing the associated uncertainties as well as the potential for these extreme outflows to affect the evolution of their host galaxies (Section 4.5).

Variation of Absorption Line Profiles
As shown in Section 3.2, each of the galaxies in our sample requires multiple components to fit the absorption troughs for the transitions studied here. The Mg II absorption troughs delineate the outflow kinematics most cleanly. Fig. 8 shows our best fit to the Mg IIλ2796 absorption trough for each galaxy and illustrates the variety of absorption profiles observed in our data.
First, there is a substantial variation in the distribution of the Mg II absorption troughs in velocity space. Two galaxies (J0905 and J1125) lack any absorption from 0 to ∼2000 km s −1 . One galaxy (J2116) shows two distinct troughs separated by ∼600 km s −1 . Five galaxies (J0826, J0901, J1232, J1558, and J2140) exhibit complex kinematics with contiguous blueshifted absorption from ∼900−1500 km s −1 . The remaining six galaxies (J0944, J1219, J1341, J1450, J1506, and J1613) have  even more complex absorption profiles with nearly continuous blueshifted absorption out to ∼2000 km s −1 . Second, there is a large variation in the number of velocity components required to describe each profile, as well as their widths (or Doppler parameter, b). As discussed in Section 3.2 we fit the blended absorption lines with the minimum number of components needed to characterize the velocity asymmetry. This number varies from two to ten velocity components. Their widths vary from 8 to 344 km s −1 , with a mean value of 106 km s −1 . Third, there is a large variation of the Mg II covering fraction (C f ) in different galaxies. For four galaxies (J0944, J1450, J1506, and J2140) most of their velocity components have C f = 1, such that the shape of their absorption profiles is determined by the optical depth rather than C f . For the rest of the sample, we find C f as low as 0.12, with a mean value of 0.57. Fourth, three galaxies (J1232, J1341, and J1450) additionally have redshifted velocity components, indicating infalling gas. Lastly, 9 of the 14 galaxies in our sample exhibit clear signs of Mg II emission with veloc-ity shifts from z sys that vary from 0 km s −1 to +450 km s −1 . Overall, there is a remarkable variation in the profiles and kinematics of these galaxies.
The galaxies in our sample are characterized by extreme and "bursty" star formation episodes that likely drive the powerful outflows observed, which can reach far into the circumgalactic medium (CGM) of the galaxy (Rupke et al. 2019). As mentioned in Section 1 our team observed the galaxy Makani with the Keck CosmicWeb Imager (KCWI; Morrissey et al. 2018) and uncovered two distinct outflows traced by [OII] emission: a largerscale, slower outflow (∼300 km s −1 ) and a smaller-scale, faster outflow (∼1500 km s −1 ). The velocities and sizes of the two wind components map exactly to two recent starburst episodes that this galaxy experienced 0.4 Gyr and 7 Myr ago, determined by its star formation history (SFH). To understand if the Mg II multiple velocity components seen here in this larger sample are also connected to the highly impulsive "burstiness" of the star formation in these galaxies we next investigate potential correlations between the Mg II absorption profiles and the SFHs of each source. Figure 9 shows the SFHs for the galaxies in our sample derived using Prospector (Johnson et al. 2021) as described in Section 2.2. The mean light-weighted age of the stellar population younger than ∼1 Gyr is reported in the bottom right of each panel. The Mg IIλ2796 absorption profile model is shown for each target in the upper right inset as presented in Fig. 8. We note that the order of the galaxies in this figure is different from the others in the paper and follows their light weighted age.
Unlike Makani which had two clear bursts and two clear outflows at different velocities, we do not find a simple correspondence between the number of bursts in the SFH and the number of Mg II absorption features in this study. Galaxies exhibiting similar Mg II absorption troughs can have a variety of SFHs and vice-versa. For example, the Mg II EWs of J0905 and J1125 are substantially different, 8.01 and 1.88Å, respectively. However, their absorption troughs are similar in several re-gards as they both show no absorption from 0 to ∼2000 km s −1 and are fit using two velocity components. Their SFHs both display a burst of star formation around 30 Myr ago, but their light-weighted ages are substantially different (17 and 102 Myr, respectively) as J0905's SFH is characterized by a younger burst ∼6 Myr ago. As another example, J1341 and J1450 have comparable Mg II EWs (10 and 10.3Å) and very broad blueshifted troughs (∼2000 km s −1 ) with extremely complex kinematics fit with a close number of components (nine and ten). Nevertheless, they have SFHs that are substantially different, with inferred light-weighted ages of 14 and 107 Myr, respectively. In particular, J1341 shows a starburst in the last 10 Myr, as opposed to J1450 which has no significant bursts in the same time frame. On the other hand, J0826 has a SFH that is similar to J1341; however J0826 exhibits a remarkably different Mg II absorption trough, with a narrower profile (∼1500 km s −1 ), no absorption at the systemic redshift (z sys ), no redshifted components, and a minimum absorption at a much higher velocity of ∼−1200 km s −1 . Additionally, J0901 has a Mg II absorption profile that is comparable to J0826 in shape and EW but has a very different SFH, with older light-weighted age (54 Myr) and a decreasing trend of star formation in the last 10 Myr. The galaxies with the fastest velocity components are not necessarily the ones with the most recent bursts, and the galaxies with strong recent bursts do not all have fast outflow components, though some do (e.g. J0826, J0905, J0944, J1219, J1341, J1613). Furthermore, we do find that eight galaxies (J0826, J0901, J0905, J1125, J1450, J1558, J1613, and J2116) have in their SFH an older burst (>10 Myr ago), regardless of the presence of a younger burst. Among these 8 galaxies, 6 (all except J0905 and J1125) have Mg II components at low velocities that may be tracing slower outflowing gas driven by the older burst of star formation, similarly to Makani.

Absorption and Emission at the Systemic Redshift
The Mg II doublet is the most sensitive transition among the absorption lines covered in this study, however, it is affected by resonantly-scattered wind emission which we expect to fill in the absorption profiles around z sys Prochaska et al. 2011). As Mg II is a resonant transition, the emitted photons are continually reabsorbed due to the absence of finestructure splitting. This trapping process interferes with the escape of the photons, complicating the interpretation of the origin of Mg II emission. In the traditional model of a galactic-scale outflow expanding as a shell, this creates a P-Cygni-like profile for each Mg II doublet component, with blueshifted absorption and redshifted emission. The Mg II absorption arises from gas moving toward the observer as it absorbs photons in its rest frame. The emission arises from the receding component of the outflow, where photons emitted in the restframe, having scattered to escape towards the observer, are redshifted and go through any intervening gas. The emission line profile is centered near the systemic redshift of the galaxy. This emission can "fill in" the absorption within ∼200 km s −1 of z sys , decreasing the EW by up to 50%, shifting the centroid of the absorption lines by tens of km s −1 , and reducing the opacity near systemic. A study of cool gas outflows that ignores this line emission may underestimate the true optical depth and/or incorrectly infer that the wind partially covers the source ).
In our sample, we detect Mg II emission in 9/14 of the galaxies, with velocity shifts of 0 km s −1 to +450 km s −1 from z sys . In J0905 we can observe emission in both Mg II lines, with an observed ratio of 1. This ratio is in the optically thick regime (e.g. Chisholm et al. 2020), and it agrees with the ratio observed in Makani (Rupke et al. 2019). The Mg II trough in J0905 is so blueshifted (∼−2400 km s −1 ) that the Mg II absorption profile is not affected by emission filling. For the remaining eight galaxies we adopt a ratio of 1 in our fits, as they clearly show Mg IIλ2803 emission but the corresponding Mg IIλ2796 line is suppressed by Mg IIλ2803 absorption at the same wavelengths. If the ratio of the two Mg II components is closer to 2 rather than 1 we may underestimate the effect of line-filling. However, the Mg II absorption profiles in these galaxies are not substantially affected by emission filling as the resonance emission is not expected to fill in the high velocity components of Mg II absorption.
In our sample the majority of the Mg II absorption EW is blueshifted; the minimum intensity of the trough lies near the systemic velocity in only two objects (J1450 and J1232). The intensity minima are blueshifted by ∼300 km s −1 to ∼2,000 km s −1 in the rest of the sample, with little or no absorption at the systemic velocity. We use the Mg IIλ2803 absorption profile to quantify the EW at systemic, as it does not suffer from blending with Mg IIλ2796. We find that 11/14 galaxies in our sample have less than four percent of the Mg II EW within 200 km s −1 of z sys . Among the objects that display Mg II emission, J1232 is the most potentially affected by emission filling as ∼26% of its Mg II EW is within 200 km s −1 of z sys . However, emission filling is not a concern for the bulk of our sample of highly blueshifted Mg II absorption lines.
The lack of substantial Mg II absorption at the systemic velocity is an interesting feature in our sample, as first noted in Perrotta et al. 2021 and Davis et al., (submitted). One interpretation of the lack of absorption at systemic could be that the extreme outflows in these galaxies may have expelled the bulk of the interstellar medium (ISM) in these sources. However, as many of the galaxies in our sample had a burst of star formation in the last 10 Myr, it is unlikely that these recent winds have entirely removed the ISM on such a short time scale. Additionally, some galaxies likely have ongoing star formation, such that the ISM can not be entirely absent.
Another possibility is that the lack of absorption at systemic is due to an observational selection effect. Our parent sample is characterized by an extremely high outflow detection rate (∼90%; Davis et al., (submitted)). While it is possible that this high outflow incidence reflects a very wide opening angle of ubiquitous outflows in these galaxies, it could also be that the magnitude and color cuts used to select our sample may have identified galaxies where a powerful outflow has excavated a hole in the ISM, causing the galaxies to appear very bright and blue. As a consequence, there may be little or no ISM left along the particular lines of sight studied here, while there is still remaining ISM in the galaxies along other sightlines.
An example of this scenario is provided by the galaxy J0905. Geach et al. (2014) use the IRAM Plateau de Bure interferometer to study the CO(2-1) emission line in J0905, which is a tracer of the bulk of the cold molecular gas reservoir. They observe a CO emission line at the systemic velocity of the galaxy, which traces ∼65% of the total cold molecular gas in the galaxy with an inferred mass of M H2 = (3.1±0.6)×10 9 M . This suggests that along some sightlines the amount of gas swept up by the outflow can be large, despite the galaxy retaining a substantial amount of its ISM.
Another potential explanation for the exceptionally high velocity Mg II absorption components and the lack of absorption at z sys is provided by strong radiative cooling (e.g., Bustard et al. 2016;. It is possible for a cold gas phase to form "in-situ" within a large-scale galactic wind via thermal instabilities and condensation of a fast-moving hot wind rather than being entrained and gradually accelerated. The byproduct of this mechanism is cold gas at similar high velocities as the hot phase. In this model, the cold clouds accelerated by the hot wind are rapidly destroyed on small scales and at low velocities by hydrodynamical instabilities. As a result, the hot wind with enhanced mass loading and density perturbations can cool radiatively on larger scales forming an extended region of atomic and ionized gas moving at ∼10 3 km s −1 , while the gas at low velocities is not observable. We come back to the origin and formation of cold gas in galactic outflows in the next Section 4.3.

Comparison with Theoretical Models
The existence of very fast, cool gas observed in outflowing winds from star-forming galaxies has been a persistent puzzle. In this Section, first, we briefly describe the processes most commonly invoked for cool gas acceleration in winds, then we use the results from this study and our parent sample to understand what insights can be gained into the role of these mechanisms in the extreme galactic outflows observed here.

Mechanisms of Cool Gas Acceleration
The cool gas phase in winds is commonly explained as the acceleration of clouds from the host ISM via ram pressure from the hot phase (e.g., Veilleux et al. 2005). However, several simulations have challenged this explanation, demonstrating that ram pressure alone is not effective at accelerating cool gas clouds to the velocities and large scales observed without the clouds being shredded by hydrodynamical instabilities and becoming incorporated into the hot flow (Cooper et al. 2009;Mc-Court et al. 2015;Scannapieco & Brüggen 2015;Schneider & Robertson 2015, 2017Zhang et al. 2017). Recent work has shown that under the right background conditions and when sufficient large clouds are considered, cool gas can survive as a result of a mixing and cooling cycle. This may increase the cool gas flux as the hot gas condenses out, effectively growing the clouds rather than destroying them (Armillotta et al. 2016;Gritton et al. 2017;Gronke & Oh 2018Fielding & Bryan 2022).
In an alternative model (Efstathiou et al. 2000;Silich et al. 2003;) the cold phase can form "in-situ" via thermal instabilities and condensation from the hot wind on large scales, provided it is sufficiently mass-loaded via the destruction of cool gas in the inner regions of the flow (see also Lochhaas et al. 2018). Additional models for cold cloud acceleration have also been proposed including momentum deposition by supernovae, the radiation pressure of starlight on dust grains (Murray et al. , 2010(Murray et al. , 2011Hopkins et al. 2012;Zhang & Thompson 2012;Krumholz & Thompson 2013;Davis et al. 2014;Thompson & Krumholz 2016), and cosmic rays (Everett et al. 2008;Socrates et al. 2008;Uhlig et al. 2012). It has also been suggested that several of these mechanisms may be taking place simultaneously in order to drive cool outflows efficiently (Hopkins et al. 2012;Veilleux et al. 2020), making it difficult to isolate the different processes potentially at play.
From an energetic point of view, starburst galaxies with powerful winds, like those in our sample, are ideal candidates for outflows driven by the radiation pressure from Eddington-limited star formation (Diamond-Stanic et al. 2012;Geach et al. 2014;Rupke et al. 2019;Perrotta et al. 2021). Recently, our team obtained results that confirm this hypothesis. Rupke et al., (submitted) present Keck/ESI (Echellette Spectrograph and Imager; Sheinis et al. 2002) long-slit spectra of the two wind episodes observed in the galaxy Makani, drawn from the same parent sample as the galaxies in this study. They infer momentum and energy outflow rates in the inner (R II = 0 − 20 kpc), recent (7 Myr ago), fast (∼2,000 km s −1 ) outflow that implies a momentumdriven flow driven by the hot ejecta and radiation pressure from the extreme, possibly Eddington-limited, compact starburst.
Here, in the Sections below, we focus on the kinematics of low-ionization absorbers tracing the cool, ionized phase of the extreme outflows in our sample to gain insights into the distribution of the outflowing gas and the physical mechanisms that may occur at the interface between the hot and cool wind phases.

Inferred Spatial Distribution of Ionized Gas
Simulations indicate that halos of Mg II emission are a ubiquitous feature across the galaxy population to at least z = 2 (Nelson et al. 2021). Mg II halos extend from a few to tens of kpc at the highest masses, i.e. far beyond the stellar component of a galaxy. Moreover, they are highly structured, clumpy, and asymmetric. DeFelippis et al. (2021) generate mock Mg II observations from the TNG50 simulation (Nelson et al. 2021) and produce absorption spectra to compare with observed data. Although the mock sight lines are too small to be comparable with observations, they show the non-uniform distribution of the Mg II absorption, usually concentrated in discrete clumps. The mock velocity spectra, despite originating from a remarkably small fraction of the total Mg II gas in the halo, reflect the diversity of the Mg II gas distribution and kinematics in halos of similar galaxy mass. On the other hand, halos with similar morphology exhibit similar mock spectra. The sample of star-forming galaxies in DeFelippis et al. (2021) differs from ours as it is characterized by galactic disks, while our galaxies are late-stage mergers with spherical morphologies. However, we also find great diversity in the Mg II absorption profiles (see Section 4.1). This suggests that the Mg II halos of our galaxies, despite spanning only 0.76 dex in stellar mass, may have diverse morphologies.
Mg II absorption traces cool ∼10 4 K gas that is usually found in simulations in high-density clouds or filaments permeating the volume-filling, low-density, hot phase at large radii (Schneider et al. 2020;Nelson et al. 2021). The remarkably extended velocity distribution of the Mg II absorption profiles in our sample may re- We find that C f determines the shape of most of the absorption troughs in 10/14 galaxies in our sample. In particular, six galaxies (J0905, J1613, J2116, J1341, J1232, and J1125) show deeper troughs at velocities where low-ionization gas covers a higher fraction of the continuum. In the remaining four galaxies (J1450, J0944, J1506, and J2140) the absorption profiles are shaped by the optical depth. Additionally, in our sample, the C f does not show a unique trend with outflow velocity.
flect the filamentary structure of dense outflowing material. Additionally, some simulated halos exhibit fountain flows in Mg II emitting gas, with signatures of infalling gas clouds in addition to wind-driven outflows (Nelson et al. 2021). We note that we find clear signatures of inflowing gas in the spectra of three galaxies in our sample (J1232, J1341, and J1450), where the components that trace infalling gas are redshifted ∼200−300 km s −1 (from the systemic redshift) and are notably narrower that the blueshifted components that trace outflowing gas. Studies of infalling gas typically utilize surveys of a hundred or more galaxies, as the detection rate of redshifted absorption lines is around 3−6% (e.g. Sato et al. 2009;Martin et al. 2012;Rubin et al. 2012). The fraction of galaxies with inflows in our sample is 21±11%, higher than in other studies. However, this value is based only on three sources, so we can not conclusively state whether the accretion rate in our sample is significantly higher.
Simulations showing the survival of cool clouds traveling through a hot medium find that the clouds undergo hydrodynamical instabilities that create an elon-gated shape with a wake (e.g., Armillotta et al. 2016;Gronke & Oh 2018). The coolest and densest gas is typically located inside the head of the cloud, while the ionized gas is found more in the turbulent wake behind the cloud, produced by the mixing between the cool gas ablated from the cloud and the hot medium. One way to investigate this potential structure observationally is to compare the absorption troughs of two transitions of the same species that have different ionization, such as neutral and singly-ionized Mg. As mentioned in Section 3.2.2, one of our approaches to fit Mg I absorption profiles is to adopt the same kinematic components as Mg II and estimate the covering fraction (C f ) that best fits the Mg I spectral region. An advantage of constraining the kinematics in this way is that we can directly compare the C f of the two different ions as a function of velocity. We find that the shallow troughs of Mg I relative to Mg II require a lower C f for the former at every velocity, with C f (Mg I) ranging from 0.16 to 0.66 C f (Mg II), with a median value of C f (Mg I) = 0.4 C f (Mg II). The similar kinematics but systematic off-set in C f implies that the Mg I absorption arises from denser regions within more extended structures traced by Mg II. This finding is in agreement with theoretical models predictions that the coolest and densest gas is typically located in the more internal and self-shielded part of the clouds.
We note that even if we do not adopt the same kinematics for Mg II and Fe II, we find a remarkable agreement in most of the velocity components. For these components, we find that C f (Fe II) is systematically lower than C f (Mg II). This is expected as Fe II traces denser gas than Mg II, and it is in line with the previous result of Mg I being less spatially extended than Mg II.

Covering Fraction Trends with Outflow Velocity
Earlier studies have tried to interpret line profiles of low-ionization ions (e.g. Mg II, Mg I, Fe II), and in particular their C f distribution, in the context of driving mechanism models for galactic winds. Martin & Bouché (2009) study a sample of five starburst galaxies at z∼0.2 and find that C f decreases as the outflow velocity increases beyond the velocity of minimum intensity, which in their case corresponds to the velocity of maximum C f . The authors interpret the velocity-dependent C f in their sample as a result of geometric dilution associated with the spherical expansion of a population of absorbers. In the context of this simple physical scenario, their result implies that the high-velocity gas detected in absorption is at a larger radii than the lower velocity (and higher C f ) gas, implying an accelerating wind. This hypothesis is also suggested by Chisholm et al. (2016) in a study of the wind-driving mechanisms and distribution of C f in a nearby starburst galaxy. This simple accelerating model does not take into account, however, the complexity of the interaction between cool clouds and the hot surrounding wind, such as radiative cooling, cloud compression due to shocks, and effects of shear flow interactions that produce hydrodynamic instabilities. More recent studies have shown that cold clouds are unlikely to survive that kind of acceleration over time (Scannapieco & Brüggen 2015;Brüggen & Scannapieco 2016;Schneider & Robertson 2017;Zhang et al. 2017). Figure 10 shows the fitted Mg II(C f ) as a function of velocity for the blueshifted (i.e. outflowing) components in our sample. C f does not show a unique trend with velocity. Seven galaxies (J0901, J0905, J1219, J1232, J1341, J1613, and J2116) have decreasing C f with increasing velocity, three galaxies (J0826, J1125, and J1558) have increasing C f with increasing velocity, and four galaxies (J0944, J1450, J1506, and J2140) have a constant C f value. This variation in the C f with outflow velocity, along with the variation of the absorption profiles described in Section 4.1, may capture the complex morphology of ∼kpc scale, inhomogeneously distributed, clumpy gas and the intricacy of the material in the turbulent mixing layers between the cold and the hot phases (e.g., Fielding et al. 2020;Nelson et al. 2021).
Our starburst galaxy sample shows that a model of an outflow accelerating over time is unlikely to be valid. One piece of evidence is the observation of the large scale, multi-phase, and multi-burst wind in the galaxy Makani. Rupke et al. (2019) find that the inner (R II = 0 − 20 kpc) and younger (∼7 Myr ago) wind is faster (with maximum speeds exceeding 2000 km s −1 ) than the more extended (R I = 20−50 kpc) and older wind (∼400 Myr ago), which has speeds of ∼100 km s −1 . Rupke et al. (submitted) find that the larger wind is consistent with having slowed down in the extended halo and CGM of the galaxy.
A second piece of evidence against an accelerating model comes from the analysis of the Mg II EW and velocity as a function of the light-weighted age of the galaxies in our parent sample. In this work, we do not find a correlation between the total Mg II EW and galaxy light-weighted age. However, the galaxies all have young light-weighted ages, between 13 and 192 Myr, and represent many of the youngest objects in the parent sample studied by our team (Davis et al., (submitted)). Using the parent HizEA sample, which has a larger dynamic range in light-weighted age, Davis et al. (in prep) find an inverse correlation between Mg II absorption EW and light-weighted age, as well as a complementary inverse correlation between outflow velocity and light-weighted age. These results are consistent with models in which the outflow decelerates with time, which support scenarios in which the cold gas condenses out of a hot flow. In particular, the observed trend of lower outflow velocity with larger light-weighted age matches the analytic models of Lochhaas et al. (2018) well, where impulsive bursts of star formation driven winds slow down and cool as they expand into the CGM.
Regardless of their physical origin, the large velocity width of the absorption troughs observed here requires contributions from multiple structures along the sightline. It is likely that our sightlines intersect many mixing layers, making it difficult to interpret the C f inferred from our fits in terms of the C f of single clouds that grow over time entrained in the hot wind (e.g., Gronke & Oh 2018;Fielding & Bryan 2022). To complicate the interpretation of the absorption line profiles further there is the multi-burst nature of the outflows in our sample. It is possible that an initial burst of star formation drives an outflow that consumes most of its energy shock-heating the surrounding ISM and consequently slows down while excavating a hole through which we observe the host galaxy as blue and bright. An outflow driven by a second burst of star formation may have a very different evolution as it can follow the path of minimum density and inertia created by the first burst, retaining more energy and velocity. A single sightline does not provide a full picture of the entire galaxy, and it can intersect multiple outflow episodes, which can overlap in the projected observed velocity space, making it challenging to distinguish the two wind episodes through absorption studies.

Comparison of MgII and FeII Absorption Troughs
It is currently challenging to use simulations to interpret observed absorption line profiles, as the models do not typically create spectra drawn from realistic sightlines through the simulated halos. They also typically identify the cool wind phase by temperature rather than making predictions for individual ions. However, we can consider their main results and model trends to make comparisons with observed spectra. Schneider et al. (2020) present a new simulation part of the Cholla Galactic OutfLow Simulations (CGOLS) project, a series of high-resolution global disk simulations of galaxy outflows. The authors model an M82-like galaxy with a 5 pc resolution in a 10 kpc volume, which is ideally suited to capture both the launching of the wind and the detailed phase structure during the subsequent expansion. Their stellar feedback injection mechanism allows the driving of energetic hot winds from a high surface density galaxy disk that gives rise to a complex, multiphase outflow. They find that the cold phase is primarily in embedded clouds that are gradually shredded and thereby enhance the mass loading of the hot phase. The authors demonstrate that the mixing between hot and cool gas in the wind is an effective way of transferring momentum from one phase to another, and this occurs at all radii. In cases where the mixed gas has a high enough density to cool, it does so with a higher velocity, leading to a linear relationship between mixed fraction and velocity.
To test this prediction we can assume that the Fe II in our spectra primarily traces the entrained cool gas phase component (i.e., gas rich in Fe results from SN Ia explosions), while the Mg II primarily traces the cooled hot wind material (i.e., gas ejected from SN II and rich in α-elements). As we do not expect the hot wind to be composed entirely of SN II ejecta, it should be α-element enriched but quite diluted compared to pure SNe ejecta. Therefore rather than searching for a complete absence of Fe II, we can instead examine how the Mg II to Fe II ratio varies as a function of outflow velocity. We expect the EW(Mg II) to EW(Fe II) ratio to increase for gas with less mixing (i.e., with less cool gas entrained as the hot wind travels through the ISM); such gas likely originates from direct cooling of the hot wind phase.
In Figure 11 we overlay the Mg IIλ2796 (blue spectrum) and Fe II λ2382 (red spectrum) absorption profiles for the galaxies in our sample. We find that the Mg IIλ2796 troughs typically have higher absorption EW. Eight galaxies (J0826, J1232, J1341, J1450, J1506, J1558, J2116, and J2140) exhibit a clear lack of Fe II at the highest velocities. To verify that this lack of Fe II is not due to a sensitivity effect, we use Equation 5 to infer the Fe II optical depth for the velocity components where Fe II is not detected. For each component, we use the optical depth of the corresponding Mg II and assume Mg II and Fe II to have the same kinematics. Equation 5 takes into account the difference in oscillator strength and dust depletion between Fe II and Mg II. As Fe II typically has lower C f compared to Mg II, we adopt a lower C f (Fe II), in line with the other Fe II absorption lines in the same spectrum. We find that in four of these galaxies (J0826, J1341, J1506, and J1558) we can not exclude that Fe II absorption may be present but too weak to detect. In the other four galaxies (J1232, J1450, J2116, and J2140) the lack of Fe II at the highest velocities is real and can be explained by gas cooling directly from the hot wind at the highest velocities. We also find that in two galaxies (J0944, J1613) there is a lack of Fe II at the lowest velocities that is not due to sensitivity effects. Finally, J2140 is the only galaxy in our sample where Fe II is not detected at any velocity. This is also the only galaxy in our sample that hosts a faint AGN (Sell et al. 2014), which could result in different feedback mechanisms and spectral features.
In Figure 12 we show the Mg IIλ2796 to Fe II λ2382 EW ratio as function of outflow velocity. To calculate this ratio we divide the Mg IIλ2796 (corrected for Mg IIλ2803 absorption) and Fe II λ2382 spectral regions in bins of 200 km s −1 and integrate the spectra, in order to estimate upper limits where the lines are not detected. The size of the circles in this figure is proportional to the strength of the EW(Mg II). Filled circles show the velocity bins where we detect Fe II, while empty circles show Fe II non-detections. This EW ratio does not have a consistent trend with the outflow velocity. However, for the eight galaxies discussed above that lack Fe II at the highest velocities, the EW ratio is larger than in the velocity bins where we detect Fe II.
In summary, several of the galaxies in our sample do not have Fe II detected at the highest outflow velocities, implying a lower mixing fraction of entrained cool gas. While the cold gas outflow phase is most likely produced  Figure 11. Comparison of the Mg IIλ2796 (blue spectrum) and Fe II λ2382 (red spectrum) line profiles in the spectra of our galaxies. The shaded blue and red regions represent the Mg IIλ2796 and Fe II λ2382 best fit profiles, respectively. The Mg IIλ2796 troughs generally have higher absorption equivalent width than the Fe II λ2382 troughs.
by a combination of physical mechanisms, these results suggest that the cold gas at the highest velocities could directly condense out of the hot wind phase, as suggested by some theoretical models.

Trends of Equivalent Width with Galaxy Properties
We now investigate the relationship between the outflow absorption strength and galaxy properties in our sample. Figure 13 compares the total Mg II equivalent width (EW) with host galaxy stellar mass (M * ), star formation rate (SFR), and star formation rate surface density (Σ SF R ). Since the Mg II doublet components are blended together in most of our galaxy spectra, we consider the total EW measured from the Mg IIλλ2796, 2803 transitions. We compare our sources with data from five representative star-forming galaxy samples in the literature that provide information about both Mg II doublet component EWs (Weiner et al. 2009;Kornei et al. 2012;Rubin et al. 2014;Bordoloi et al. 2014;Prusinski et al. 2021).
To quantify the relations presented in Figure 13, we use the python Markov Chain Monte Carlo package PyMC3 (Salvatier et al. 2016) to compute a linear fit and the associated uncertainties. In each panel, we show the linear fits with 1 and 3σ error regions shaded in blue. We also characterize the correlations using the Spearman rank correlation coefficient (ρ) and the Pearson linear correlation coefficient (R). These values and their corresponding statistical significance are listed in the figure caption.
A strong positive correlation between EW(Mg II) and M * is evident in panel (a) of Figure 13, while EW(Mg II) exhibits a weaker correlation with SFR in panel (b) and Σ SF R in panel (c), as has been found in previous work (e.g., Martin et al. 2012;Rubin et al. 2014). Our galaxies follow the general trends seen in other samples, though with substantial scatter. The EW(Mg II) versus M * correlation has been explained in these previous works as being due to two factors: 1) the EW of ISM absorption (at the systemic redshift) increases with increasing M * , and 2) emission filling is stronger in lower mass galaxies, preferentially suppressing absorption around the systemic redshift. As a result, the lower emission filling in the spectra of massive galaxies (relative to lower mass galaxies) amplifies the positive correlation between EW(Mg II) and M * . However, for the eight galaxies (J0826, J1232, J1341, J1450, J1506, J1558, J2116, and J2140) that have a clear lack of Fe II at the highest velocities, the EW ratio is larger than in the velocity bins where we detect Fe II, suggesting that we are seeing gas directly condensed out from the hot wind.
As described in Section 4.2, our sample of massive galaxies lacks substantial ISM absorption at the systemic velocity, and the large EWs(Mg II) observed here are typically from absorption that is blueshifted more than 200 km s −1 , such that emission filling is not substantial. The systematically larger EWs of massive galaxies may be due to a larger number of blueshifted components compared to lower mass galaxies, which would imply a correlation between the maximum velocity of the outflows and M * . This has been observed (e.g., Davis et al., submitted, Rubin et al. 2014;Heckman & Borthakur 2016), though the correlation with stellar mass does not appear to be as strong as with other galaxy properties. Therefore, the explanation previously proposed for the EW(Mg II) vs. M * correlation is plausible.
Another potential contributing factor to the EW(Mg II) vs. M * correlation is that as M * increases, the reservoir of gas that can potentially be part of the outflow should increase as well, for star-forming galaxies. In this context, the mass and velocity of the outflowing gas are consequences of the star formation conditions in the galaxy. There is growing agreement that Σ SF R is one of the most important properties governing the velocities of galactic winds (e.g., Davis et al., submitted, Heckman et al. 2015;Heckman & Borthakur 2016;Prusinski et al. 2021). If we consider two galaxies with similar EW(Mg II), M * , and SFR but different Σ SF R , we may expect the galaxy with higher Σ SF R to have a higher outflow velocity, as the similar energy released by the supernovae explosion is injected into a smaller volume creating a more explosive event. We see this trend if we consider the most massive starburst galaxies from Rubin et al. (2014) and Kornei et al. (2012) with EW(Mg II), M * , and SFR comparable to our sample. The ∼2 dex larger Σ SF R in our sample is reflected in the substantially higher outflow velocities observed (Davis et al., submitted). A specific example is J1232, which has Σ SF R similar to the most massive galaxies from Rubin et al. (2014) and Kornei et al. (2012) and it also has a comparable outflow velocity v 98 = 761 km s −1 (v avg = 217 km s −1 ), which is the lowest in our sample. In , and ΣSF R (c) for the galaxies presented in this work (filled circles; see legend in the upper left corner in the panel (a)) compared with star-forming galaxies with clear wind signatures: 0.3 < z < 1.4 galaxies from Rubin et al. (2014) shown with blue stars; z ∼ 1 galaxies from Kornei et al. (2012) shown with blue circles; 1 ≤ z ≤ 1.5 galaxies from Prusinski et al. (2021) shown with blue rectangles; stacked spectra of 1,406 galaxies at z ∼ 1.3 − 1.5 from Weiner et al. (2009) shown with dark blue diamonds; coadded spectra of 486 galaxies at 1 < z < 1.5 from Bordoloi et al. (2014) shown with dark blue reversed triangles. Error bars are included only for results using stacked data, for clarity. Fit lines are an MCMC-generated linear fit to the data, with 1 and 3σ error regions shaded in blue. The Spearman rank order correlation coefficient ( this context, considering two galaxies with similar M * , SFR, and Σ SF R but very different EW(Mg II), we may expect that the galaxy with lower EW(Mg II) would have a higher outflow velocity, as a comparable amount of energy injected has to accelerate less gas. An example of this is J1125, which has similar M * , SFR, and Σ SF R to J1450, but has ∼5 times lower EW(Mg II) and a substantially larger outflow velocity of v 98 = 2244 km s −1 (v avg = 1813 km s −1 ), compared to v 98 = 1874 km s −1 (v avg = 529 km s −1 ) for J1450. While the general trends described above apply to many galaxies in Figure 13, there are exceptions, and these relations have substantial scatter. While such scaling relations are not trivial to interpret, they can be useful from a statistical standpoint. However, when we look at the details of individual galaxies, such general trends may not strictly apply, especially given that single sightlines can provide only a partial view of the entire galaxy. Another complication is that EWs measured from saturated absorption lines reflect only a lower limit on the wind absorption strength.

Mass Outflow Rates
In this section, we estimate mass outflow rates for the observed winds in these galaxies and discuss the assumptions used to determine these values. We note that the mass outflow rates are confined to the phase we probe with our observations, which is the ionized gas phase, and therefore they are lower limits on the total mass outflow rates which include neutral and molecular gas (e.g. Bolatto et al. 2013;Chisholm et al. 2016).
Following Rupke et al. (2005b) we assume a simple model for the wind that depends on the physical parameters derived from our line fitting. If we consider a single thin shell wind at radius r, with thickness d r, the average mass outflow rate is: where µm p is the mean atomic weight (m p is the proton mass and µ = 1.4 is the correction for the relative He abundance). We calculate time-averaged outflow rates because they are a more useful quantity than instantaneous values to use when comparing to star formation rates. To average the mass outflow rate over the wind lifetime and obtain Equation 6, we divide the instanta-neousṀ by t wind = r/v i . The sum in Equation 6 is performed over the individual outflowing (i.e., blueshifted) velocity components in each galaxy: Ω i is the solid angle subtended by a given component as seen from the wind's origin, N i (H) is the total hydrogen column density of that component, and v i is the central velocity of that component. In this model, we split Ω (the wind's global covering factor) into two parts to account for a potential biconical morphology and a clumpy wind, rather than assuming a smooth shell. We model Ω in terms of the large-scale covering factor C Ω , related to the wind's opening angle, and the local covering factor C f , related to the wind's clumpiness. Thus, Ω/4π = C Ω C f . We adopt C Ω = 1 based on the high outflow detection rate in the parent sample (Davis et al., (submitted)). We can then rewrite Equation 6 aṡ An estimate of the total hydrogen column density in the outflow requires knowledge of the ionization state and metallicity of the gas, as well as the amount of dust depletion for the element employed to derive N (H). Our spectral coverage provides access to a series of strong Fe II resonance lines which have oscillator strengths spanning a substantial range and therefore can be used to estimate the column density of singly ionized iron. From this one can determine the total hydrogen column density for an assumed metallicity and ionization fraction. However, as we detect Mg II in a larger number of absorption components for most galaxies in our sample, we begin with estimates of the total hydrogen column derived from Mg II. We describe below how these values might need to be adjusted to provide more accurate estimates using results from Fe II where they exist.
We adopt χ(Mg II) = 0.7 since, in the case of Mg, the singly-ionized state Mg II is likely dominant in photoionized gas at T ∼10 4 K (Murray et al. 2007;Churchill et al. 2003). We assume a solar abundance ratio (log Mg/H = −4.45; Asplund et al. 2021), which is consistent with results from an ensemble of line ratio diagnostic diagrams used to estimate the metallicities of galaxies in our sample (Perrotta et al. 2021). We assume a dust depletion factor (d(Mg)) of −0.5 dex for Mg, as measured in the local Galactic ISM (Jenkins 2009). The total hydrogen column density in the outflow is where the sum is performed over the individual Mg II outflowing velocity components in each galaxy. We list in Table 2 the N (H) values derived using this method. These values are lower limits not only because of our assumed conservative ionization correction, dust depletion, and abundance ratios but also because the absorption in the Mg II components is saturated for most of the galaxies in our sample. We emphasize that most of the saturated Mg II absorption troughs in our spectra do not appear black which implies that C f < 1. Because the actual hydrogen column densities could be substantially higher than those estimated above from Mg II, we can use bounds on τ 0 (Mg II) from Equation 5 to obtain a revised N (Mg II) for all the Mg II components that are also detected in Fe II. Using these values, we derive updated N (H) that we report in Table 2. In our galaxy sample the hydrogen column density for the same component obtained from Equation 8 using τ 0 (Mg II) meas (i.e. measured directly from Mg II fits) and τ 0 (Mg II) inf (i.e. inferred using bounds from Fe II fits) differ from a factor ranging from 0.3 to 19. The total hydrogen column density for each galaxy is obtained as the sum of N (H) for the individual velocity components. Using the highest N (H) values in Table 2 (i.e., obtained using τ 0 (MgII), or τ 0 (FeII)), we obtain total hydrogen outflowing column densities for our galaxies of N (H) tot = 3 ×10 19 − 2 ×10 21 cm −2 , with a median column density of 2 ×10 20 cm −2 . For most galaxies in our sample, these values are still lower limits, as Fe II provides a bound on τ 0 only for 36% of the Mg II absorption troughs in our sample. Moreover, our estimates are derived assuming a conservative depletion of iron (relative to magnesium) onto dust grains. Dust depletion factors (d(X)) measured in the local Galactic ISM fall in the range −(1.0−2.3) dex for Fe (Jenkins 2009). Here we adopt d(Fe) = −1 dex; a greater depletion correction would increase τ 0 (Mg II) meas , and consequently N (H), by a factor of ∼3 − 60. There is therefore an order of magnitude uncertainty due to the dust depletion correction alone.
Returning to our goal of estimating the mass outflow rates in the winds observed in absorption, we can use the values constrained by our fits for all of the variables in Equation 7 except for the spatial extent of the wind (r). While the absorption lines are sensitive to absorbing gas at any location along the line of sight to the observed galaxies, they do not provide information on the outflow geometry. Here we choose to adopt a thin shell of uniform radius 5 kpc. This radius is motivated both by observations of star-forming galaxies at similar redshifts (e.g., Burchett et al. 2021;Zabl et al. 2021) and by our own integral field spectrograph data of outflows observed in our sources (Rupke et al. 2019 and Mg II nebulae reaching far beyond the stars in the galaxies, with radial extents of a few tens of kpc (Perrotta et al. in prep). Therefore, a physical extent of 5 kpc for the outflows in this paper is used as a reasonable order of magnitude estimate. Determining a more accurate value is challenging as these galaxies show variations object to object and ion to ion. SinceṀ is directly proportional to r, any change in the physical extent of the outflow will proportionally result in a change toṀ . For example, if the value of r increases by a factor of two,Ṁ will increase by a factor of two. Assuming a thick wind instead of a thin wind will decrease the derived outflow rates since the radial factor in Equation 7 is the inner radius in the thick wind case, rather than the outer radius. For example, a thick wind extending from 1 to 5 kpc has a mass 5 times lower than the r = 5 kpc thin shell that we adopt here. We prefer a shell geometry over a constant velocity wind given the morphology of the ionized outflow observed in Makani that resembles an evacuated and limb-brightened bipolar bubble. Moreover, a constant velocity wind is likely not appropriate for our sample, which has had recent strong star formation bursts due to major merger activity. Finally, the absorption spectra presented here do not support a constant wind velocity model given their multi-component nature and large velocity extent. We emphasize that while the geometry of the outflow (e.g., physical extent, thickness) is uncertain, this is subdominant to other uncertainties such as the dust depletion factor and optical depth, as discussed above. We may rewrite Equation 7 using fiducial values as follows: TheṀ values for our sample derived using Equation 7 are reported in Table 5. TheṀ values obtained using constraints only from Mg II fits are listed in column 1, while column 2 listsṀ calculated utilizing constraints from Mg II and, when available, Fe II fits. We have robust estimates (i.e., non lower-limits) on the ionized gas column density for only two sources in our sample, J0905 and J1219. For these two objects using the constraints from the Fe II measurements increases the mass outflow rates values from their lower limits to values that are roughly an order of magnitude larger. This underscores that the lower limits derived for the rest of the galaxies in our sample may underestimate the actual mass outflow rates by a factor of 10.
As described above, deriving accurateṀ estimates requires precise knowledge of the physical conditions in the outflowing gas, the properties of the interstellar medium where the outflow propagates, and the outflow geometry and kinematics. We discussed above the assumptions adopted here to determine these first-order estimates, and we discussed how they represent lower limits on the actual outflow rates for most of the galaxies in our sample. We stress that the biggest uncertainty in calculating  Figure 14. Estimated ionized gas phase mass outflow rates versus star formation rates for our galaxies. These estimates represent lower limits on the actual ionized mass outflow rates for the whole sample, except for J0905 and J1219. For these two galaxies using the constraints from the Fe II measurements increases the mass outflow rates values from their lower limits to values that are roughly an order of magnitude larger. The tilted lines (from bottom to top) illustrate mass outflow rates that are 0.1, 1, and 5 times the star formation rate.
the mass outflow rates is introduced by the uncertainties on the N (H) estimates. In our study, the leading uncertainty on N (H) arises from the difficulty of breaking the degeneracy between τ and C f for saturated Mg II doublet lines, when additional constraints from the Fe II multiplet are not available. We also report in Table 5 the mass loading factor, η ≡Ṁ /SFR. η compares the amount of gas entrained by the outflow to the amount of gas actively being converted into stars, and it represents a useful way to quantify how the outflow may impact the future evolution of the host galaxy, as it relates to the rate at which gas is being removed from the galaxy to the rate at which it is currently creating stars. Figure 14 shows the mass outflow rates as calculated using Mg II only constraints (Table 5, column 1) as a function of the star formation rate. As discussed these are lower limits on the total ionized gas mass outflow rates. In order to demonstrate the potential magnitude of the systematic underestimate ofṀ derived using Mg II only, we show theṀ values calculated using bounds on τ 0 from Fe II fits (Table 5, column 2) for the two sources (J0905 and J1219) in our sample that have the most robust estimates on τ 0 and which are therefore not lower limits. The twoṀ estimates for J0905 and J1219 are connected by lines, to facilitate comparison. There is a large correction for these two sources between the Mg II only estimates (lower limits) and the more robust estimates (not lower limits), of factors of ∼15 and 12, respectively. The rest of the galaxies in the sample might have a similar magnitude correction if better constraints on N (H) from the Fe II absorption lines were available. Applying similar corrections to the rest of the sample (i.e. to theṀ MgII values reported in Table 5) as those obtained for J0905 and J1219, we would infer mass outflow rates in the range ∼50 − 2000 M yr −1 for our galaxies. Figure 14 also includes lines of constant mass loading factors, η. The two galaxies with more robust constraints onṀ , J0905 and J1219, exhibit η values of the order of ∼25 and 17, respectively. Such high values of the mass loading factor indicate that the out-flows are ejecting gas at a much higher rate than the production rate of stars. Similar η values are found in ultraluminous infrared galaxies (ULIRGs) at z < 0.5 (Rupke et al. 2005c). While the mass outflow rates are first-order estimates only, they are informative to verify that the observed starburst-driven outflows can clearly affect the evolution of their host galaxies by actively removing substantial amounts of gas. The mass outflow rates for J0905 and J1219 are estimated to be ∼2200 and 1500 M yr −1 , which are exceptionally high mass outflow rates. This finding is consistent with results based on observations of molecular gas for two galaxies (J0905 and J1341) in our sample (Geach et al. 2014(Geach et al. , 2018. In particular, Geach et al. (2014) use the IRAM Plateau de Bure interferometer to study the CO(2-1) emission line, a tracer of the bulk of the cold molecular gas reservoir, in the galaxy J0905. They find that a third of the total molecular gas content appears to have been ejected on a scale of approximately 10 kpc. They estimate that the kinetic energy associated with the outflowing component is consistent with the momentum flux available from stellar radiation pressure, demonstrating that nuclear starbursts are capable of ejecting large amounts of cold gas from the central regions of galaxies.
As discussed in Section 4.3.4, the cold gas outflow phase is most likely produced by a combination of multiple physical processes. Therefore the high mass loading factors derived here do not necessarily refer to a single production mechanism of the cool gas. However, as a reference, we compare our estimates with predictions from Lochhaas et al. (2021), who model hot supersonic winds driven by supernovae energy injection in star-forming galaxies. The authors derive the characteristic momentum rates of hot galactic winds that can undergo single-phase cooling on large scales and produce cold gas as a function of the SFR of the galaxy. Interestingly, the predicted velocities of the wind at the maximum wind momentum rate for SFRs similar to our galaxies are compatible with the highest velocity absorption components in our spectra. However, most of our sample shows higher mass loading factors than the theoretical maximum predicted by their model. This would imply that our observations do not trace only cool outflows that obtained their momentum directly from the free-flowing hot wind, but rather that there may be some mass loading outside the wind-driving region. This supports the hypothesis that the cold gas outflow phase may have multiple origins.
In summary, when robust constraints from the Fe II multiplet are available we find extreme mass outflow rates (∼1500 − 2200 M yr −1 ) and η values of ∼17 − 25. If similar corrections are applied to the rest of the galax-ies, they would haveṀ in the range ∼50−2000 M yr −1 and ∼90% of the sample would have η > 1.

SUMMARY AND CONCLUSIONS
We use new optical Keck/HIRES spectroscopy of 14 compact starburst galaxies at z ∼0.5 to probe the small-scale structure of the powerful galactic outflows observed in our sample, gain insights into the role of various physical mechanisms in these extreme feedback episodes, and investigate their potential impact on the evolution of their host galaxies. These galaxies are massive (M * ∼10 11 M ), extremely compact (half-light radius ∼few hundred pc), have very high star formation rates (mean SFR ∼ 200 M yr −1 ) and star formation surface densities (mean Σ SFR ∼ 2000 M yr −1 kpc −2 ), and host extremely fast (mean maximum velocity ∼−2000 km s −1 ) outflows traced by Mg II absorption lines (Tremonti et al. 2007;Davis et al., (submitted). The high resolution ( 8 km s −1 ) spectra presented here cover a suite of Mg and Fe absorption lines (Mg IIλλ2796, 2803doublet, Mg Iλ2852, and Fe IIλ2344, 2374, 2582, 2587, 2600. This exquisite data set allows us to directly measure precise column densities and covering fractions (C f ) as a function of velocity for these ions and characterize the kinematics of the cool gas outflow phase (T ∼10 4 K). Our main conclusions are as follows: 1) Mg II absorption troughs best delineate the outflow kinematics among the ions studied here. We find a substantial variation in the absorption profiles in our data (Section 4.1 and Figure 8). In particular, there is a large variation in the minimum number of velocity components required to fit each profile, ranging from two to ten; in their widths, Doppler parameter b, ranging from 8 to 344 km s −1 with a mean of 106 km s −1 ; and in their C f , with four galaxies having C f = 1, while the rest of the sample has C f as low as 0.12 and a mean value of 0.57.
2) We investigate if the multiple velocity components observed in our sample are related to the extreme and "bursty" star formation episodes in these galaxies. There is not a simple correlation between the number of bursts in the star formation history (SFH) and the Mg II absorption profiles (Section 4.1 and Figure 9). Galaxies showing similar Mg II absorption troughs can have a variety of SFHs and vice-versa.
3) Mg II emission is detected in 9/14 of the galaxies in our sample with velocity shifts of 0 km s −1 to +450 km s −1 from z sys . Emission filling is not an issue for most of our galaxies, which have highly blueshifted Mg II absorption lines (Section 4.2 and Figure 8). Indeed, 11/14 galaxies have less than four percent of the EW(Mg II) within 200 km s −1 of z sys , where emission line filling can be important. We also find a lack of substantial Mg II absorption at the systemic velocity, and we present three possible explanations for this: the bulk of the interstellar medium (ISM) is expelled by the powerful outflows, biased observational selection criteria, and/or Mg II traces a cold gas phase formed "insitu" within a galactic wind via thermal instabilities and condensation of the fast-moving hot phase. 4) Mg II, Mg I, and Fe II exhibit remarkably similar absorption profiles, suggesting these species reside in the same, low-ionization gas structures. However, Mg II has on average a higher C f at a given outflow velocity than Fe II, and in all cases a higher C f than neutral Mg, implying that the absorbing clouds are not homogeneous (Section 4.3.2). This result is in agreement with theoretical predictions that the coolest and densest gas is typically located in the more internal and self-shielded part of the clouds. 5) We find that C f does not display a unique trend with velocity (Section 4.3.3 and Figure 10). This variation, along with the variation of the absorption profiles, may capture the complex morphology of ∼kpc scale, inhomogeneously-distributed, clumpy gas, and the intricacy of the material in the turbulent mixing layers between the cold and the hot phases. Moreover, other observations of our galaxies are consistent with models in which the outflows decelerate with time as they expand into the CGM. 6) We consider the possibility of Fe II as the primary tracer of the entrained cool gas phase component and Mg II as the tracer of the cooled hot wind material and examine how their ratio varies as a function of outflow velocity. Several of the galaxies in our sample do not have Fe II detected at the highest outflow velocities, indicating a lower mixing fraction of entrained cool gas. This suggests that the cold gas at the highest velocities in these galaxies most likely directly condensed out of the hot wind phase, as suggested by some theoretical models (Section 4.3.4, Figure 11 and 12). 7) We estimate mass outflow rates for the observed winds in our sample and discuss how these estimates are lower limits on the actual outflow rates for most of the galaxies (Section 4.5 and Figure 14). We show that the two galaxies in our sample that have robust constraints from the Fe II multiplet have extremely high mass outflow rates (∼ 1500 − 2200 M yr −1 ) and mass loading factors (η ∼17−25). If similar corrections are applied to the rest of the galaxies we inferṀ ∼50−2200 M yr −1 with a typical value of η ∼5. This demonstrates that starburst galaxies are capable of ejecting very large amounts of cold gas which will substantially impact their future evolution.
The galaxy sample studied here provides a prime opportunity to investigate star formation and feedback at its most extreme. In a forthcoming paper based on integral field unit Keck/KCWI spectra, we will focus on studying the morphology, physical extent, and kinematics of the outflows in our sample. Such data on these galaxies provides strong observational constraints to theoretical simulations that aim to implement realistic stellar-driven galactic outflows. The comparison of outflow properties between simulations and observations will advance our understanding of galactic feedback, especially from stellar processes, during a critical stage of massive galaxy evolution. from the Heising-Simons Foundation grant 2019-1659. S. P. and A. L. C. acknowledge support from the Ingrid and Joseph W. Hibben endowed chair at UC San Diego. The data presented herein were obtained at the W. M. Keck Observatory, which is operated as a scientific partnership among the California Institute of Technology, the University of California, and the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation. The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.

A. NOTES AND FITS OF INDIVIDUAL TARGETS
Figures 15−24 show the normalized spectra of the other 10 galaxies in our sample centered on the transitions relevant to this study: Mg IIλλ2796,2803Å, Mg Iλ2803Å, Fe IIλ2382Å, Fe IIλ2600Å, Fe IIλ2344Å, Fe IIλ2586Å, and Fe IIλ2374Å. For each galaxy, we show the Mg Iλ2853 best fit profile as derived following the constrained approach. Details for each galaxy are given below. Figure 15 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J0826+4305 (hereafter J0826). J0826 spectrum shows Mg II emission within −400 and +400 km s −1 of z sys . In this galaxy, we see strong absorption from Mg II falling within −1560 and −75 km s −1 of z sys . The Mg II trough is well described by eight velocity components that show complex kinematics. The Mg IIλ2796 and Mg IIλ2803 troughs exhibit comparable intensity at all velocities for seven over eight components. In that case, the Mg II traces optically thick gas. As the lines are not black the covering fraction (C f < 1) determines the shape of the absorption troughs. The most blueshifted component (v = −1404 km s −1 ) has τ 0 < 1 and C f = 1. In this galaxy, we note that the C f increases with increasing blueshift from z sys . We detect Fe II and Mg I absorption within −1350 and −1000 km s −1 of z sys , which we model using two velocity components. They have quite good kinematic correspondence (within the errors) to two Mg II components. However, the Fe II and Mg I (independent) fit find the least blueshifted component to be substantially broader than the corresponding Mg II (20 and 55%, respectively). We find both Fe II components to trace optically thin gas and have lower C f than Mg II (32 and 45% of C f (Mg II), respectively). Based on comparison of τ 0 (Mg II) meas and τ 0 (Mg II) inf , we conclude that the Mg II fit gives a lower limit estimate of the N(Mg II) for the two components we have detection of both Mg II and Fe II. The actual τ 0 (Mg II) is ∼3 times larger than the measured one. As for previous targets we favor the Mg I constrained fit solution because the neutral Mg fraction must be less than one percent to produce an optically thin Mg Iλ2853 trough with C f =1 as found by the independent fit. The constrained fit finds for both Mg I components a smaller C f than Mg II (28 and 33% of C f (Mg II), respectively). Figure 16 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J0901+0314 (hereafter J0901). J0901 spectrum exhibits Mg II redshifted emission (+93 km s −1 ) within −300 and +500 km s −1 of z sys . In this galaxy, we detect strong absorption from Mg II falling within −1750 and −100 km s −1 of z sys . The Mg II trough is well characterized by four velocity components that exhibit complex kinematics. One Mg IIλ2796 component blends together with three Mg IIλ2803 ones. Our Mg II best fit includes one optically thin component (at v = −591 km s −1 ) with C f =1, and three more blueshifted components. The latters trace optically thick gas, with Mg IIλ2796 and Mg IIλ2803 showing nearly identical intensity at all velocities. Their absorption profiles are not black and their shape is determined by C f < 1. We detect Fe II and Mg I absorption within −1550 and −1000 km s −1 of z sys , which we model using two velocity components. We find that one of these two components has a good kinematics correspondence (comparable v and b to within the error bars) to one Mg II component. C f for both Fe II and Mg I (from the constrained fit), is found to be lower than C f (Mg II), 52% and 43% of C f (Mg II), respectively. The values of τ 0 (Mg II) meas and τ 0 (Mg II) inf for the component detected at −1266 km s −1 represent a lower limit (the actual τ 0 (Mg II) is ∼3 times larger), while agreeing well within the errors for the component detected at −1426 km s −1 . The independent Mg I fit shows kinematics in line with the constrained fit. However, to be a valid solution, the neutral Mg fraction must be less than one percent. The second and more blueshifted Fe II and Mg I (from the independent fit) component in our model show a substantially broader absorption profile (46% and 72% of b(Mg II), respectively). Based on comparison of τ 0 (Mg II) meas and τ 0 (Mg II) inf we infer that Mg II should be optically thin, but that is inconsistent with the nearly identical intensity of Mg IIλ2796 and Mg IIλ2803. We attribute this discrepancy to the low SNR in the Fe II and Mg I profiles blueward of −1350 km s −1 that is inadequate to definitively constrain the kinematics of the absorption troughs. Figure 17 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J0944+0930 (hereafter J0944). J0944 spectrum does not show Mg II emission. In this galaxy, we see strong absorption from Mg II falling within −2100 and −120 km s −1 of z sys . The Mg II trough is well described by five velocity components that show convoluted kinematics. Most of them (four Mg IIλ2796, and four Mg IIλ2803) blend together. We find that the profile shape for all Mg II components is determined by the optical depth, rather than the C f , that is unity. Our Mg II best fit includes only one optically thick component (at v = −1125 km s −1 ), the remaining components trace optically thin gas. The Mg I trough has a remarkable alignment with Mg II. We model the Mg I absorption using five components, as well. In the constrained fit the five Mg I components are described by lower C f compared to Mg II, with values in the range of 32 − 56% of C f (Mg II). The independent Mg I fit finds a great agreement in terms of kinematics with the corresponding Mg II components. However, to be a valid solution, the neutral Mg fraction must be less than one percent. Therefore, we favor the constrained fit result. We detect Fe II within −1350 and −700 km s −1 of z sys , that we model using three velocity components. They have good alignment with three Mg II components, but they are characterized by broader line profiles (in the range of 30 to 50 %). Based on a comparison of τ 0 (Mg II) meas and τ 0 (Mg II) inf we conclude that the Mg II fit provides a decent estimate of the N(Mg II). However, we note that the difference in the kinematics of the two absorption troughs makes the comparison of the τ 0 not accurate.

J1125−0145
Figure 18 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J1125−0145 (hereafter J1125). J1125 spectrum does not exhibit Mg II emission. In this galaxy, we see absorption from Mg II falling within −2170 and −1830 km s −1 of z sys . The Mg II trough is well described by two optically thick velocity components. The identical intensity at all velocities of the Mg II absorption lines indicates that their profile shapes are determined by C f < 1. We detect Mg I and Fe II within the same velocity range as Mg II. We model the Mg I and Fe II troughs using two components and find a remarkable kinematic correspondence (comparable v and b to within the errors) to the Mg II fit. Both, Mg I (from constrained fit) and Fe II absorption profiles are described well by lower C f than Mg II(∼50 and ∼43 % of C f (Mg II), respectively). Comparing τ 0 (Mg II) meas and τ 0 (Mg II) inf we argue that the Mg II fit provides a lower limit estimate of the N(Mg II) and that the actual central optical depth for the two Mg IIλ2803 components is closer to 9, 22 (rather than 2 and 3 as measured from the Mg II fit). Figure 19 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J1219+0336 (hereafter J1219). J1219 spectrum very clearly shows Mg IIλ2803 emission within −600 and +600 km s −1 of z sys . The corresponding Mg IIλ2796 line is not obvious because of Mg IIλ2803 absorption at the same wavelengths. The inclusion of an emission component in the model improves the overall fit to the Mg II absorption trough. In this galaxy, we see strong absorption from Mg II falling within −2400 and −150 km s −1 of z sys . The Mg II trough is well described by four velocity components that blend together producing extremely intricate kinematics. Our Mg II best fit includes one optically thin component (at v = −755km s −1 ) with C f =1. The remaining three components trace optically thick gas, and their profiles are shaped by C f < 1. We detect Mg I within −2400 and −950 km s −1 of z sys . We model the Mg I trough using three components that show a good qualitative alignment with Mg II. The constrained Mg I fit finds lower C f than Mg II, in the range 29 − 56 % of C f (Mg II). In our independent Mg I fit model, we find for the most blueshifted component an excellent kinematic agreement to the corresponding Mg II one (v and b within the errors). The other two components show kinematics substantially different from Mg II. We attribute this difference to the low SNR in the Mg I profile redward of −1400 km s −1 which is not sufficient to conclusively constrain the kinematics of the absorption troughs. We detect Fe II within −2400 and −750 km s −1 of z sys . Before performing the Fe II fit, we identify and model the Mn II λλλ 2576, 2594, and 2606 triplet (yellow solid lines). We model the Mg I trough using three components that show extremely good kinematic agreement with Mg II. We find Fe II to have lower C f than Mg II (53, 46, and 55 % of C f (Mg II), respectively). Based on a comparison of τ 0 (Mg II) meas and τ 0 (Mg II) inf , we conclude that the Mg II fit gives only a lower limit to the N(Mg II), and that the actual τ 0 for the two Mg IIλ2803 components, is closer to 90, 21, and 91 (than 10, 1.2, and 10). Figure 19 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J1506+5402 (hereafter J1506). J1506 spectrum very clearly shows redshifted (∼126 km s −1 ) Mg IIλ2803 emission from −300 to +550 km s −1 of z sys . The corresponding Mg IIλ2796 line is not obvious because of Mg IIλ2803 absorption at the same wavelengths. The inclusion of an emission component in the model improves the overall fit to the Mg II absorption trough. In this galaxy, we see strong absorption from Mg II falling within −2030 and +150 km s −1 of z sys . The Mg II trough is well characterized by five velocity components that blend together creating complex kinematics. We find that all Mg II components trace optically thin gas, and their absorption profiles are shaped by the optical depth, rather than the C f , that is unity. We detect Mg I within −1750 and −750 km s −1 of z sys . We model the Mg I trough using two components that show a good alignment with Mg II. The constrained fit finds lower C f (Mg I) than Mg II (45 and 60% of C f (Mg II), respectively). Adopting the independent Mg I fit we find the two Mg I components to be slightly narrower (7 and 20%) than the corresponding Mg II. In both cases, the Mg I traces optically thin gas as Mg II. However, we favor the constrained fit solution since to have C f (Mg I) = 1 the neutral Mg fraction must be less than one percent. We detect Fe II within −1650 and −750 km s −1 of z sys . We model the Fe II troughs using two components that show a good alignment with Mg II. However, we note that b(Fe II) is 38% narrower than b(Mg II) for the most blueshifted component, and 15% broader for the other one. Despite the discrepancy between their kinematics, we find a very good agreement between τ 0 (Mg II) meas and τ 0 (Mg II) inf within the errors. Figure 19 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J1558+3957 (hereafter J1558). J1558 spectrum distinctly shows slightly redshifted (57 km s −1 ) Mg IIλ2803 emission within −650 and +700 km s −1 of z sys . The corresponding Mg IIλ2796 line is not apparent because of Mg IIλ2803 absorption at the same wavelengths. The inclusion of an emission component in the model improves the fit to the Mg II absorption trough. In this galaxy, we see strong absorption from Mg II falling within −1350 and −100 km s −1 of z sys . The Mg II trough is well described by four velocity components that exhibit intricate kinematics. Our Mg II best fit finds the most blueshifted component to be optically thing and with C f = 1. For the other three components, the measured equivalent width ratio for Mg IIλ2796 and Mg IIλ2803 is not consistent (∼30 to 65% lower) with the optically thin limit. They trace optically thick gas and their absorption profiles are shaped by C f < 1. We detect weak Mg I within −1350 and −630 km s −1 of z sys . We model the Mg I trough using two components. The constrained fit finds Mg I to be characterized by lower C f than Mg II (31 and 50% of C f (Mg II), respectively). The independent Mg I fit finds two components that show a good qualitative alignment with Mg II, but fairly broader absorption profiles (56 and 18% respectively). We prefer the constrained fit solution since to have optically thin Mg I for the least blueshifted component and C f (Mg I) = 1 the neutral Mg fraction must be less than one percent. We detect Fe II within −1000 and −630 km s −1 of z sys . We model Fe II using only one component with a good kinematic correspondence to one Mg II component. We find that Fe II traces optically thin gas and the profile is well described by a lower C f than Mg II(51 % of C f (Mg II)). Based on comparison of τ 0 (Mg II) meas and τ 0 (Mg II) inf , we conclude that the Mg II fit gives a lower limit estimate of the N(Mg II) for the component we have detection of both Mg II and Fe II. The actual τ 0 (Mg II) is ∼4 times larger than the measured one. Figure 19 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J1613+2834 (hereafter J1613). J1613 spectrum clearly exhibits redshifted (284 km s −1 ) Mg IIλ2803 emission within −700 and +1200 km s −1 of z sys . The corresponding Mg IIλ2796 line is not obvious because of Mg IIλ2803 absorption at the same wavelengths. The inclusion of an emission component in the model improves the fit to the Mg II absorption trough. In this galaxy, we see strong absorption from Mg II falling within −2570 and −50 km s −1 of z sys . The Mg II trough is well described by five velocity components that blend together and produce complicated kinematics. We find that Mg II traces optically thick gas. The measured equivalent width ratio for Mg IIλ2796 and Mg IIλ2803 is not consistent (∼25 to 45% lower) with the optically thin limit. The absorption profile of all five components is well described by C f < 1. We do not identify any significant absorption trough in the Mg I spectral region. We detect weak Fe II in correspondence with the most blueshifted Mg II component. We model the Fe II trough using only one component. Despite the Fe II component showing a good alignment with Mg II (similar v within the errors), our best fit finds a significantly (44%) broader absorption profile. We attribute this to the low SNR in the Fe II spectral region. According to our best fit Fe II traces optically thin gas. We infer from the Fe II fit that the corresponding Mg II component should be less optically thick, with τ 0 (Mg II) inf = 2.7 and τ 0 (Mg II) meas = 10. This is probably due to a difference in the kinematics (b(Mg II) = 57% b(Fe II)) of the two absorption lines that makes the comparison of the τ 0 not accurate, underestimating τ 0 (Mg II) inf .

J2116−0624
Figure 23 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J2116−0624 (hereafter J2116). J2116 spectrum does not exhibit Mg II emission. In this galaxy, we identify two Mg II absorption troughs at around −280 and −1430 km s −1 of z sys , respectively. We find that the most blueshifted Mg II component (v = −1428 km s −1 ) traces optically thick gas and its profile shape is characterized by C f <1. The least blueshifted Mg II component (v = −284 km s −1 ) is close to the transition between optically thin and thick gas, with τ 0 =0.97. The measured equivalent width ratio for Mg IIλ2796 and Mg IIλ2803 is ∼25% lower than the theoretical value. Therefore, we conclude that the absorption profile shape is determined by the C f <1 rather than the optical depth. We detect weak Mg I and Fe II in correspondence with the least blueshifted Mg II component. We model Mg I and Fe II using only one component. The independent Mg I fit finds a substantially broader (∼60%) absorption profile than the corresponding Mg II. We attribute this discrepancy to the difficulties of identifying such a weak absorption trough. As for previous targets we favor the Mg I constrained fit solution because the neutral Mg fraction must be less than one percent to produce an optically thin Mg Iλ2853 trough with C f =1 as found by the independent fit. The constrained fit finds Mg I to be characterized by a smaller C f than Mg II (32% of C f (Mg II)). The Fe II fit shows a very good kinematic agreement with Mg II (similar v and bwithin the errors). The Fe II absorption traces optically thin gas with C f =1. We find a good agreement between τ 0 (Mg II) meas and τ 0 (Mg II) inf within the errors. Figure 24 shows our best fit to the Mg II, Mg I, and Fe II absorption lines in the galaxy J2140+1209 (hereafter J2140). J2140 spectrum does not exhibit Mg II emission. In this galaxy, we identify strong Mg II absorption troughs falling within −950 and +250 km s −1 of z sys . The Mg II trough is well described by five velocity components that marginally blend together (one Mg IIλ2796 and three Mg IIλ2803 components). We find that Mg II traces optically thick gas for the least blueshifted component (v = −48 km s −1 ). The measured equivalent width ratio for Mg IIλ2796 and Mg IIλ2803 is not consistent (∼40% lower) with the optically thin limit. The absorption profile of this component is shaped by C f < 1. The remaining four components trace optically thin gas and their profile shape is well described by the optical depth. We do not identify any significant absorption trough in the Mg I and Fe II spectral regions.  Figure 16. J0901 spectrum, which exhibits Mg II emission within −300 and +500 km s −1 of zsys. In this galaxy, we see strong absorption from Mg II falling within −1750 and −100 km s −1 of zsys. The Mg II trough is well described by four velocity components. We detect Fe II and Mg I absorption within −1550 and −1000 km s −1 of zsys, which we model using two velocity components.  Figure 17. J0944 spectrum, which does not exhibit Mg II emission. In this galaxy, we see strong absorption from Mg II falling within −2100 and −120 km s −1 of zsys. The Mg II trough is well described by five velocity components. The detected Mg I trough has a remarkable alignment with Mg II. We model the Mg I absorption using five components. We detect Fe II within −1350 and −700 km s −1 of zsys, which we model using three velocity components.  Figure 19. J1219 spectrum, which shows Mg II emission within −600 and +600 km s −1 of zsys. In this galaxy, we see strong absorption from Mg II falling within −2400 and −150 km s −1 of zsys. The Mg II trough is well described by four velocity components. We detect Mg I within −2400 and −950 km s −1 of zsys. We model the Mg I trough using three components that show a good qualitative alignment with Mg II. We detect Fe II within −2400 and −750 km s −1 of zsys. We model the Fe II trough using three components that show good kinematic agreement with Mg II.  Figure 20. J1506 spectrum, which shows Mg II emission within −300 and +550 km s −1 of zsys. In this galaxy, we see strong absorption from Mg II falling within −2030 and +150 km s −1 of zsys. The Mg II trough is well described by five velocity components. We detect Mg I within −1750 and −750 km s −1 of zsys. We detect Fe II within −1650 and −750 km s −1 of zsys. We model the Mg I and Fe II troughs using two components that show a good qualitative alignment with Mg II. . J1558 spectrum, which shows Mg II emission within −650 and +700 km s −1 of zsys. In this galaxy, we see strong absorption from Mg II falling within −1350 and −100 km s −1 of zsys. The Mg II trough is well described by four velocity components. We detect Mg I within −1350 and −630 km s −1 of zsys. We model the Mg I trough using two components that show a good alignment with Mg II. We detect Fe II within −1000 and −630 km s −1 of zsys. We model Fe II using one component.  Figure 22. J1613 spectrum, which shows Mg II emission within −700 and +1200 km s −1 of zsys. In this galaxy, we see strong absorption from Mg II within −2570 and −50 km s −1 of zsys. The Mg II trough is well described by five velocity components. We do not detect Mg I. We detect Fe II at about −2400 km s −1 of zsys. . J2140 spectrum, which does not exhibit Mg II emission. In this galaxy, we see strong absorption from Mg II falling within −950 and +250 km s −1 of zsys. The Mg II trough is well described by five velocity components that marginally blend together. We do not detect any Mg I and Fe II absorption lines.