# The Circumgalactic Medium of Extreme Emission Line Galaxies at
$\boldsymbol{z\sim 2}$: Resolved Spectroscopy and Radiative Transfer Modeling
of Spatially Extended Ly$\boldsymbol{\alpha}$ Emission in the KBSS-KCWI
Survey111Based on data 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 NASA, and was made possible by
the generous financial support of the W.M. Keck Foundation.

Dawn K. Erb The Leonard E. Parker Center for Gravitation, Cosmology and
Astrophysics, Department of Physics,
University of Wisconsin-Milwaukee, 3135 N Maryland Avenue, Milwaukee, WI
53211, USA  Zhihui Li (李智慧) Cahill Center for Astronomy and Astrophysics,
California Institute of Technology, 1200 E California Blvd, MC 249-17,
Pasadena, CA 91125, USA  Charles C. Steidel Cahill Center for Astronomy and
Astrophysics, California Institute of Technology, 1200 E California Blvd, MC
249-17, Pasadena, CA 91125, USA  Yuguang Chen (陈昱光) Cahill Center for
Astronomy and Astrophysics, California Institute of Technology, 1200 E
California Blvd, MC 249-17, Pasadena, CA 91125, USA  Department of Physics and
Astronomy, University of California Davis, 1 Shields Avenue, Davis, CA 95616,
USA  Max Gronke Max-Planck Institute for Astrophysics, Karl-Schwarzschild-Str.
1, D-85741 Garching, Germany  Allison L. Strom Department of Physics and
Astronomy and Center for Interdisciplinary Exploration and Research in
Astrophysics (CIERA), Northwestern University, 2145 Sheridan Road, Evanston,
IL 60208, USA  Ryan F. Trainor Department of Physics and Astronomy, Franklin &
Marshall College, 637 College Avenue, Lancaster, PA 17603, USA  Gwen C. Rudie
The Observatories of the Carnegie Institution for Science, 813 Santa Barbara
Street, Pasadena, CA 91101, USA

###### Abstract

The resonantly scattered $\rm Ly\alpha$ line illuminates the extended halos of
neutral hydrogen in the circumgalactic medium of galaxies. We present integral
field Keck Cosmic Web Imager observations of double-peaked, spatially extended
$\rm Ly\alpha$ emission in 12 relatively low-mass ($M_{\star}\sim 10^{9}$ $\rm
M_{\odot}$) $z\sim 2$ galaxies characterized by extreme nebular emission
lines. Using individual spaxels and small bins as well as radially binned
profiles of larger regions, we find that for most objects in the sample the
$\rm Ly\alpha$ blue-to-red peak ratio increases, the peak separation
decreases, and the fraction of flux emerging at line center increases with
radius. We use new radiative transfer simulations to model each galaxy with a
clumpy, multiphase outflow with radially varying outflow velocity, and self-
consistently apply the same velocity model to the low ionization interstellar
absorption lines. These models reproduce the trends of peak ratio, peak
separation and trough depth with radius, and broadly reconcile outflow
velocities inferred from $\rm Ly\alpha$ and absorption lines. The galaxies in
our sample are well-described by a model in which neutral, outflowing clumps
are embedded in a hotter, more highly ionized inter-clump medium (ICM), whose
residual neutral content produces absorption at the systemic redshift. The
peak ratio, peak separation and trough flux fraction are primarily governed by
the line-of-sight component of the outflow velocity, the H I column density,
and the residual neutral density in the ICM respectively. Azimuthal
asymmetries in the line profile further suggest non-radial gas motions at
large radii and variations in the H I column density in the outer halos.

Galaxy evolution (594), High-redshift galaxies (734), Circumgalactic medium
(1879), Galaxy spectroscopy (2171), Radiative transfer simulations (1967)

††journal: AAS Journals††software: KCWI
DRP222222https://github.com/kcwidev/kderp, astropy (Astropy Collaboration et
al. 2013, 2018), dynesty (Skilling 2004, 2006; Speagle 2020), spectral-cube
(Ginsburg et al. 2019), SExtractor (Bertin & Arnouts 1996), cmasher (van der
Velden 2020), seaborn (Waskom et al. 2020)

##  1 Introduction

A star-forming galaxy in the early universe is the nexus of a complex
interchange of gas between stars and a nested series of gaseous reservoirs:
the interstellar medium (ISM), consisting of the gas among the stars, the
circumgalactic medium (CGM), and finally the intergalactic medium (IGM). As
the transition region between the stars and the IGM, most of the key processes
of galaxy evolution are modulated through the CGM (see Tumlinson et al. 2017
for a recent review). Outflows powered by star formation drive gas out of the
galaxy and into the CGM, where it may be reaccreted onto the galaxy or
continue on to leave the galaxy entirely (e.g. Veilleux et al. 2020). New fuel
for star formation is accreted through the CGM, likely via dense, cold streams
of gas (Kereš et al. 2005; Dekel & Birnboim 2006). Ionizing photons that make
their way out of the galaxy must also traverse the neutral hydrogen in the CGM
(Rudie et al. 2013; Steidel et al. 2018).

The strongest emission line arising from gas in the CGM is due to the $\rm
Ly\alpha$ transition of hydrogen, and deep observations have now revealed that
the diffuse, distant universe is aglow with $\rm Ly\alpha$ emission (Wisotzki
et al. 2018; Ouchi et al. 2020). This emission arises primarily from faint
halos extending to tens of kpc around galaxies, but is also seen in the form
of larger nebulae (“blobs”, e.g. Fynbo et al. 1999; Steidel et al. 2000) and
filaments (e.g. Cantalupo et al. 2014; Umehata et al. 2019; Daddi et al.
2021). The initial detections of spatially extended $\rm Ly\alpha$ emission
surrounding typical star-forming galaxies at high redshifts came from stacked,
narrowband images (Steidel et al. 2011; Momose et al. 2014, 2016; Xue et al.
2017), but more recently the Multi-Unit Spectroscopic Explorer (MUSE, Bacon et
al. 2010) at the ESO-VLT and the Keck Cosmic Web Imager (KCWI, Martin et al.
2010; Morrissey et al. 2012) have enabled the study of individual halos around
galaxies from $2\lesssim z\lesssim 6$ (Wisotzki et al. 2016; Leclercq et al.
2017; Wisotzki et al. 2018; Erb et al. 2018; Leclercq et al. 2020; Chen et al.
2021).

A number of different mechanisms have been proposed to account for this
extended $\rm Ly\alpha$ emission. Perhaps the most straightforward of these is
the resonant scattering of $\rm Ly\alpha$ photons produced in galaxies by
neutral hydrogen in the CGM (Zheng et al. 2011; Kusakabe et al. 2019; Byrohl
et al. 2021), with the $\rm Ly\alpha$ profile then reflecting the kinematics
and geometry of the CGM gas. Other possible sources of the $\rm Ly\alpha$ halo
emission include in situ photoionization (fluorescence), either by ionizing
radiation escaping from the galaxy or by an external radiation field
(Kollmeier et al. 2010; Cantalupo et al. 2012; Mas-Ribas & Dijkstra 2016);
cooling radiation from infalling gas (Haiman et al. 2000; Dijkstra & Loeb
2009; Faucher-Giguère et al. 2010; Lake et al. 2015); or $\rm Ly\alpha$
emission from faint satellite galaxies (Mas-Ribas et al. 2017). Multiple
mechanisms may contribute in a given halo, with their relative importance
varying with radius (Mitchell et al. 2021). Recent theoretical results suggest
that the $\rm Ly\alpha$ properties of gaseous halos are primarily influenced
by galactic outflows within $\sim 50$ kpc, while cold accretion flows dominate
at larger radii (Chung et al. 2019); this result is in agreement with
observations that find a transition between outflow and inflow-dominated
kinematics at similar radius (Chen et al. 2020).

While most studies of $\rm Ly\alpha$ halos to date have focused on the spatial
distribution of the emission via imaging (either from narrowband filters or
reconstructed from integral field unit [IFU] data cubes), a number of recent
IFU studies have analyzed spectral variations in the extended emission, using
small samples of gravitationally lensed (Patrício et al. 2016; Claeyssens et
al. 2019; Solimano et al. 2022) and unlensed (Erb et al. 2018; Leclercq et al.
2020) galaxies at $z>2$. The inclusion of spectroscopic information has the
potential to be a powerful discriminant among the proposed emission
mechanisms, although the resonant nature of $\rm Ly\alpha$ emission has made
the extraction of physical quantities from observed spectra difficult, even in
the case of single, spatially integrated line profiles. Due to multiple
scatterings, the emergent profile depends on the kinematics, geometry and
density of neutral hydrogen and on the dust content (see Dijkstra 2014a for a
review). The strongest peak of the observed $\rm Ly\alpha$ profile is almost
always redshifted relative to the systemic redshift of the galaxy due to
backscattering from a receding galactic outflow, and when the opacity to $\rm
Ly\alpha$ photons in the outflow is relatively low (usually seen in lower
mass, highly ionized galaxies) a secondary, blueshifted peak may be visible as
well. In the local universe the separation between the two peaks has been
observed to correlate with the escape of ionizing Lyman continuum radiation,
with objects with narrower peak separations having higher escape fractions
(Verhamme et al. 2017; Izotov et al. 2018).

Spatially resolved spectroscopic studies of individual $\rm Ly\alpha$ halos
have so far mostly been based on MUSE data, and have therefore necessarily
focused on galaxies at $z>3$. These MUSE studies have analyzed the spectral
properties of $\rm Ly\alpha$ emitters with a single peak, finding that the
velocity shift of the line is generally smaller for higher surface brightness
regions and that there is a correlation between the width and velocity shift
of the line, with broader emission often tending to come from the outer halo
(Claeyssens et al. 2019; Leclercq et al. 2020; Solimano et al. 2022). At
$z=2.3$, Erb et al. (2018) studied a single low-mass galaxy with KCWI,
measuring variations in the peak ratio and separation of the double-peaked
$\rm Ly\alpha$ profile across the extended halo and finding that higher blue-
to-red peak ratios and narrower separations tended to be found at larger
radii. These spectroscopic studies have been broadly interpreted in the
context of the resonant scattering of $\rm Ly\alpha$ photons in a galactic
outflow, but definitive models for the observed trends have yet to be
constructed.

A number of radiative transfer (RT) codes have successfully reproduced the
$\rm Ly\alpha$ profiles of large numbers of spatially integrated spectra,
generally by modeling the outflow as a spherical, expanding shell (e.g.
Verhamme et al. 2008, 2015; Hashimoto et al. 2015; Yang et al. 2017; Gronke
2017). These models have provided constraints on the properties of the
scattering medium, while also indicating that even within the simplified
regime of the shell model the interpretation of the $\rm Ly\alpha$ profile is
complex. In general, the separation between the two peaks of the line
increases with increasing H I column density, while the relative strength of
the blue peak decreases with increasing velocity of the shell (e.g. Verhamme
et al. 2015). However, the physical parameters inferred from shell models do
not always match constraints on the gas obtained from interstellar absorption
and nebular emission lines (e.g. Kulas et al. 2012; Leitherer et al. 2013;
Orlitová et al. 2018), with the models predicting lower outflow velocities and
higher intrinsic line widths. More generally, the outflowing gas in real
galaxies is multiphase and spans a wide range in velocity, in contrast to the
single value assumed by the shell models (e.g. Steidel et al. 2010).

Alternatively, $\rm Ly\alpha$ RT has been studied in a more realistic
multiphase, clumpy medium, where cool, H I clumps are embedded in a hot,
highly ionized inter-clump medium (ICM) (e.g. Neufeld 1991; Hansen & Oh 2006;
Dijkstra & Kramer 2012; Laursen et al. 2013; Duval et al. 2014; Gronke &
Dijkstra 2016). In this multiphase, clumpy model, the kinematics, covering
factor, and column density of the clumps, along with the residual H I number
density in the ICM, act together to shape the morphology of the $\rm Ly\alpha$
profile. Such a clumpy model converges to the monolithic shell model in the
limit of being “very clumpy” (i.e. having $\sim$ 1000 clumps on average per
line of sight), but its unique flexibility offers the possibility of obtaining
more physically reasonable parameters of the gaseous medium that are
consistent with other observations (Li & Gronke 2022).

These models were first applied to fitting the KCWI-observed $\rm Ly\alpha$
profiles of several regions in the $z=3.1$ $\rm Ly\alpha$ blob SSA22-LAB1
(Steidel et al. 2000) by Li et al. (2021), who managed to reproduce the
diverse morphologies of the observed profiles with reasonable physical
parameters of the gaseous medium. Notably, they found that many of the
observed $\rm Ly\alpha$ profiles have significant residual fluxes at the line
center, which correspond to relatively few clumps per line-of-sight and low
residual H I density in the ICM. In addition, the very broad $\rm Ly\alpha$
wings can be reproduced by large random velocity dispersions of the clumps,
but are hard to explain in the context of shell models without requiring
unphysically large widths of the intrinsic profiles of the $\rm Ly\alpha$
emission.

Follow-up work by Li et al. (2022) modeled the $\rm Ly\alpha$ profiles of
another $z=3.1$ $\rm Ly\alpha$ blob, SSA22-LAB2, with both the multiphase,
clumpy models and shell models. They identified a significant correlation
between the shell expansion velocity and the clump outflow velocity, and found
that the multiphase, clumpy model may alleviate the inconsistencies between
the shell model parameters and the observational data. Moreover, for the first
time, they attempted to use radially-binned models to fit the spatially
resolved $\rm Ly\alpha$ profiles. They found that the $\rm Ly\alpha$ profiles
at different impact parameters can be reproduced self-consistently assuming a
common central source, and that the variation of the clump outflow velocity
with respect to impact parameter can be explained by a line-of-sight
projection effect of a radial outflow. In this paper, we build on the
methodology of Li et al. (2022) and continue to model spatially resolved $\rm
Ly\alpha$ spectra with the multiphase, clumpy model.

We analyze the spectral properties of spatially extended $\rm Ly\alpha$
emission for a sample of 12 relatively low-mass, low-metallicity galaxies at
$z\sim 2$, using integral field spectroscopy from KCWI. Our focus on low-mass
galaxies with extreme nebular line emission is motivated by the likely
importance of faint galaxies to reionization (e.g. Kuhlen & Faucher-Giguère
2012; Robertson et al. 2015) and by the observed and expected connections
between the $\rm Ly\alpha$ profile and Lyman continuum escape (Dijkstra et al.
2016; Verhamme et al. 2017; Izotov et al. 2018; Steidel et al. 2018). Escaping
ionizing radiation must travel through the CGM, and spatially resolved models
of extended $\rm Ly\alpha$ emission offer the possibility of obtaining
constraints on the physical conditions in the multiphase CGM gas. The double-
peaked nature of $\rm Ly\alpha$ emission from highly ionized sources also
provides additional constraints on the models; all 12 of our targets have
double-peaked profiles, which we quantify in both individual spaxels and
binned regions before modeling the results with state-of-the-art radiative
transfer codes.

We describe our sample selection, observations, and data reduction in Section
2, and measure the global properties of the $\rm Ly\alpha$ emission in Section
3. In Section 4 we quantify the $\rm Ly\alpha$ profiles across the extended
halos, measuring the line morphology in individual spaxels and small spatial
bins. We bin the data with larger regions in Section 5, to measure both
average properties and maximum and minimum gradients in peak ratio and
separation. In Section 6 we apply new models to the both the spatially
resolved $\rm Ly\alpha$ emission and the rest-frame UV interstellar absorption
lines, and we summarize our results and discuss their implications in Section
7. We assume the Planck Collaboration et al. (2020) values of the cosmological
parameters, $H_{0}=67.7$ km s-1 Mpc-1, $\Omega_{\rm m}=0.31$, and
$\Omega_{\Lambda}=0.69$; with these values, 1 arcsec subtends a distance of
8.4 proper kpc at $z=2.3$, the median redshift of our sample.

Table 1: Targets Observed

ID | RA | Dec | $\mathcal{R}$ | $M_{\rm UV}$ | $\beta$ | $z_{\rm neb}$ | $\log(M_{\star}/M_{\odot})$ | SFR | sSFR | $t_{\rm exp}$
---|---|---|---|---|---|---|---|---|---|---
| (J2000) | (J2000) | (AB mag) | (AB mag) |  |  |  | ($\rm M_{\odot}$ yr-1) | (Gyr-1) | (hr)
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11)
Q0142-BX165 | 01:45:16.867 |  $-$09:46:03.47 | 23.51 |  $-$21.62 |  $-$1.90 | 2.3577 |  19.13 | 24.4 | 18.1 | 5.0
Q0142-BX186 | 01:45:17.484 |  $-$09:45:07.99 | 25.32 |  $-$19.81 |  $-$1.24 | 2.3569 |  18.59 | 12.2 | 31.3 | 5.0
Q0207-BX87 | 02:09:44.234 |  $-$00:04:13.51 | 23.84 |  $-$21.15 |  $-$1.72 | 2.1924 |  19.48 |  18.3 |  12.8 | 4.7
Q0207-BX144 | 02:09:49.209 |  $-$00:05:31.67 | 23.75 |  $-$21.22 |  $-$2.03 | 2.1682 |  19.22 | 18.9 | 11.4 | 4.5
Q0449-BX110 | 04:52:17.201 |  $-$16:39:40.64 | 23.94 |  $-$21.17 |  $-$1.72 | 2.3355 |  19.29 | 18.2 |  19.3 | 5.0
Q0449-BX115 | 04:52:17.861 |  $-$16:39:45.36 | 24.88 |  $-$20.23 |  $-$2.28 | 2.3348 |  18.90 |  12.1 |  12.6 | 5.0
Q0821-MD36 | 08:21:11.410 |  $+$31:08:29.44 | 24.48 |  $-$20.82 |  $-$1.62 | 2.5830 |  19.12 | 24.1 | 18.3 | 5.1
Q1549-BX102 | 15:51:55.982 |  $+$19:12:44.20 | 24.32 |  $-$20.67 |  $-$1.64 | 2.1934 |  19.64 |  16.0 |  11.4 | 5.0
Q1700-BX729 | 17:01:27.773 |  $+$64:12:29.48 | 24.02 |  $-$21.14 |  $-$1.87 | 2.3993 | 10.10 | 24.5 |  11.9 | 4.3
Q2206-BX151 | 22:08:48.674 |  $-$19:42:25.42 | 23.91 |  $-$21.09 |  $-$2.10 | 2.1974 |  19.97 |  15.5 |  10.6 | 4.9
Q2343-BX418 | 23:46:18.571 |  $+$12:47:47.36 | 23.99 |  $-$21.10 |  $-$2.05 | 2.3054 |  18.68 | 14.4 | 30.0 | 4.8
Q2343-BX660 | 23:46:29.433 |  $+$12:49:45.55 | 24.17 |  $-$20.81 |  $-$1.87 | 2.1742 |  18.73 | 13.9 | 25.8 | 5.0

††footnotetext: Columns: (1) Galaxy ID; (2) Right ascension in hours, minutes
and seconds; (3) Declination in degrees, minutes and seconds; (4) Observed
$\mathcal{R}$-band AB magnitude; (5) Absolute UV magnitude at $\sim 2100$ Å;
(6) Rest-frame UV slope $\beta$ measured from $G-\mathcal{R}$ color; (7)
Systemic redshift from rest-frame optical nebular emission lines; (8) Stellar
mass from SED fit; (9) Star formation rate from $\rm H\alpha$ luminosity (see
Section 2.1); (10) Specific star formation rate SFR/$M_{\star}$; (11) Total
KCWI integration time.

Figure 1: Continuum-subtracted $\rm Ly\alpha$ profiles from spatially
integrated spectra, with double asymmetric Gaussian fits (discussed in Section
4) shown in orange.

Figure 2: Continuum-subtracted $\rm Ly\alpha$ images. White contours show the
$\rm Ly\alpha$ surface brightness, with the same levels in each panel:
$1\times 10^{-18}$ (dotted), $5\times 10^{-18}$ (dashed), $1\times 10^{-17}$
(dash-dot), and $2\times 10^{-17}$ (solid)
$\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{arcsec}^{-2}$, and
black contours indicate the adjacent UV continuum measured in a rest-frame 75
Å window redward of the $\rm Ly\alpha$ emission line.

##  2 Sample, Observations and Data Reduction

Table 2: Nebular Emission Line Measurements

ID | $F_{\rm H\alpha}$ | $F_{\rm H\beta}$ |  $\rm H\alpha$/$\rm H\beta$ | [N II]/$\rm H\alpha$ | [O III]/$\rm H\beta$ | O32 | O2 | $n_{e}$
---|---|---|---|---|---|---|---|---
| ($10^{-16}$ $\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}$) | ($10^{-16}$ $\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}$) |  |  |  |  |  | (cm-3)
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9)
Q0142-BX165 | 1.51 $\pm$ 0.102 | 0.42 $\pm$ 0.003 | 3.56 $\pm$ 0.24 | 0.025 $\pm$ 0.006 | 6.30 $\pm$ 0.06 | 3.29 $\pm$ 0.76 | 1.32 $\pm$ 0.12 | 94 $\pm$ 96
Q0142-BX186 | 0.60 $\pm$ 0.031 | 0.15 $\pm$ 0.005 | 4.04 $\pm$ 0.25 |  $<$0.054 | 7.48 $\pm$ 0.29 | 3.75 $\pm$ 0.82 | 1.31 $\pm$ 0.18 | 100 $\pm$ 144
Q0207-BX87 | 0.96 $\pm$ 0.081 | 0.37 $\pm$ 0.007 | 2.62 $\pm$ 0.23 | 0.041 $\pm$ 0.009 | 6.29 $\pm$ 0.16 | 6.45 $\pm$ 1.93 | 0.95 $\pm$ 0.12 | 572 $\pm$ 260
Q0207-BX144 | 1.53 $\pm$ 0.032 | 0.44 $\pm$ 0.024 | 3.44 $\pm$ 0.20 |  $<$0.031 | 6.15 $\pm$ 0.37 | 5.90 $\pm$ 1.16 | 1.24 $\pm$ 0.09 | 157 $\pm$ 84
Q0449-BX110 | 1.10 $\pm$ 0.014 | 0.30 $\pm$ 0.013 | 3.65 $\pm$ 0.17 |  $<$0.036 | 7.21 $\pm$ 0.33 | 3.86 $\pm$ 0.61 | 0.99 $\pm$ 0.05 | 488 $\pm$ 94
Q0449-BX115 | … | 0.13 $\pm$ 0.021 | … | … | 5.45 $\pm$ 0.92 | … | … | …
Q0821-MD36 | 1.14 $\pm$ 0.048 | 0.31 $\pm$ 0.027 | 3.66 $\pm$ 0.35 |  $<$0.063 | 6.94 $\pm$ 0.64 | 4.42 $\pm$ 1.41 | 0.89 $\pm$ 0.04 | 708 $\pm$ 112
Q1549-BX102 | 0.69 $\pm$ 0.132 | 0.29 $\pm$ 0.011 | 2.42 $\pm$ 0.47 |  $<$0.090 | 5.24 $\pm$ 0.22 | 3.87 $\pm$ 2.71 | 1.01 $\pm$ 0.13 | 448 $\pm$ 223
Q1700-BX729 | 1.57 $\pm$ 0.019 | 0.46 $\pm$ 0.026 | 3.41 $\pm$ 0.20 | 0.059 $\pm$ 0.007 | 5.78 $\pm$ 0.34 | … | … | …
Q2206-BX151aa$\rm H\alpha$ flux measurement from Keck-NIRSPEC, Kulas et al. (2012). | 0.63 $\pm$ 0.030 | … | … | … | … | … | … | …
Q2343-BX418 | 1.21 $\pm$ 0.027 | 0.39 $\pm$ 0.016 | 3.09 $\pm$ 0.14 | 0.043 $\pm$ 0.008 | 7.32 $\pm$ 0.30 | 7.26 $\pm$ 1.12 | 0.93 $\pm$ 0.07 | 609 $\pm$ 156
Q2343-BX660 | 1.63 $\pm$ 0.060 | 0.79 $\pm$ 0.032 | 2.06 $\pm$ 0.11 | 0.017 $\pm$ 0.005 | 6.93 $\pm$ 0.28 | 14.69 $\pm$ 2.77 | 1.15 $\pm$ 0.11 | 258 $\pm$ 128

††footnotetext: Columns: (1) Galaxy ID; (2) Observed $\rm H\alpha$ flux from
Keck-MOSFIRE; (3) Observed $\rm H\beta$ flux; (4) $\rm H\alpha$/$\rm H\beta$
flux ratio; (5) [N II]$\lambda$6585/$\rm H\alpha$ flux ratio, with 3$\sigma$
upper limits for [N II] non-detections; (6) [O III]$\lambda$5008/$\rm H\beta$
flux ratio; (7) O32 $\equiv$ [O III]$\lambda\lambda$4960,5008/[O
II]$\lambda\lambda$3727,3729, corrected for extinction; (8) Observed [O II]
flux ratio, O2 $\equiv$ [O II]$\lambda$3729/[O II]$\lambda$3727; (9) Electron
density from [O II] ratio.

###  2.1 Sample selection and properties

Target selection for this study was motivated by the simultaneous goals of
characterizing the CGM in relatively low-mass, extreme emission line galaxies
and improving our ability to extract physical information from double-peaked
$\rm Ly\alpha$ emission. Because low-mass, low-metallicity, and highly ionized
galaxies tend to exhibit strong, double-peaked $\rm Ly\alpha$ emission (e.g.
Henry et al. 2015; Erb et al. 2016; Trainor et al. 2016; Verhamme et al. 2017;
Matthee et al. 2021), these objectives largely lead to the same targets.

Our targets are drawn from the Keck Baryonic Structure Survey (KBSS; Rudie et
al. 2012; Steidel et al. 2014; Strom et al. 2017) of star-forming galaxies at
$z\sim 2$. Selection is primarily based on nebular emission line measurements
from KBSS-MOSFIRE (Steidel et al. 2014; Strom et al. 2017): five of our 12
targets are drawn from the sample of Erb et al. (2016), who studied the $\rm
Ly\alpha$ properties of $z\sim 2$ galaxies with extreme nebular emission line
ratios placing them in the upper left corner of the [N II]/$\rm H\alpha$ vs.
[O III]/$\rm H\beta$ “BPT” diagnostic diagram (Baldwin et al. 1981), with
$\log(\mbox{[{N~{}II}]}/\rm H\alpha)\leq-1.1$ and $\log(\mbox{[{O~{}III}]}/\rm
H\beta)\geq 0.75$. Galaxies in this region of the diagram lie at the low
metallicity, high ionization end of the star-forming sequence, and the $z\sim
2$ galaxies in our sample have typical metallicities $12+\log(\rm O/H)\approx
8.0$ (see Erb et al. 2016 for discussion).

Four additional targets (Q0142-BX186, Q0449-BX110, Q0821-MD36, and
Q1700-BX729) meet the nebular line ratio criteria used for the Erb et al.
(2016) paper but were identified later, and the remaining three objects,
Q0449-BX115, Q1549-BX102 and Q2206-BX151, were selected from among the
strongest $\rm Ly\alpha$-emitters (LAEs) in the $z\sim 2$ KBSS sample.
Q1549-BX102 lies just outside the emission line selection region, and
Q0449-BX115 and Q2206-BX151 cannot be placed on the diagram due to
insufficient data (Q0449-BX115 is detected in [O III] but not $\rm H\alpha$,
and Q2206-BX151 has only $\rm H\alpha$ observations). The median redshift of
the sample is $z_{\rm med}=2.32$, and all of the targets fall above the
canonical LAE threshold, with rest-frame equivalent width $W_{\rm Ly\alpha}$
$>20$ Å measured from long-slit spectroscopy.

The sample galaxies also have very high equivalent width [O III] $\lambda
5008$ emission, with $W_{[\rm O\,\sc{III}]}=870$ Å measured from a composite
$H$-band spectrum (we do not measure individual equivalent widths because the
continuum is noisy in many of the individual spectra); this value is
comparable to that of $z\sim 1$–2 reionization-era analogs selected for
extreme [O III] emission (Tang et al. 2019).

Global properties of the galaxies are given in Table 1, and nebular emission
line measurements in Table 2. The sample is largely blue and bright, with 75%
of the objects brighter than $M^{*}_{\rm UV}=-20.70$ at $z\sim 2.3$ (Reddy &
Steidel 2009) and median UV slope $\beta_{\rm med}=-1.87$. The median stellar
mass of the sample is $1.5\times 10^{9}$ $\rm M_{\odot}$, from modeling the
spectral energy distributions with the BPASSv2.2 stellar population synthesis
models (Stanway & Eldridge 2018) and assuming the SMC extinction law (Gordon
et al. 2003) and an initial mass function with slope $-2.35$ over the range
0.5–100 $\rm M_{\odot}$ and $-1.35$ between 0.1–0.5 $\rm M_{\odot}$. This
median stellar mass implies a halo mass of $\sim$$3\times 10^{11}$ $\rm
M_{\odot}$ (Girelli et al. 2020), roughly three times lower than the typical
halo mass of the KBSS parent sample (Adelberger et al. 2005; Conroy et al.
2008; Trainor & Steidel 2012); the corresponding virial radius is $\sim 60$
kpc. The SED modeling also indicates that the galaxies are young and
relatively unreddened, with median age 100 Myr and median $E(B-V)_{\rm
cont}=0.05$.

We use the $\rm H\alpha$/$\rm H\beta$ ratio and the SMC extinction law to
correct the nebular emission lines for internal reddening, finding median
$E(B-V)_{\rm neb}=0.14$. Two galaxies in the sample (Q0207-BX87 and
Q2343-BX660) have $\rm H\alpha$/$\rm H\beta$ less than the theoretical value,
here assumed to be 2.79 corresponding to Case B recombination at an electron
temperature of 15,000 K. These galaxies are assigned $E(B-V)_{\rm neb}=0$. For
the two objects that do not have measurements of the Balmer decrement we
instead use reddening measurements from the SED fitting, $E(B-V)_{\rm
cont}=0.0$ for both Q0449-BX115 and Q2206-BX151.

Star formation rates are computed from the dust-corrected $\rm H\alpha$
luminosity using the calibration of Theios et al. (2019), who calculate the
conversion between SFR and $\rm H\alpha$ luminosity for the BPASSv2.2 stellar
models used for the SED fitting described above. The resulting SFRs range from
2 to 25 $\rm M_{\odot}$ yr-1, with a median of 14 $\rm M_{\odot}$ yr-1. The
galaxies with the two lowest SFRs in the sample, Q0449-BX115 and Q2206-BX151,
are also the two for which we have determined the reddening using results from
the SED fitting; $E(B-V)_{\rm cont}$ is typically smaller than $E(B-V)_{\rm
neb}$, so for these two objects we have potentially underestimated the
extinction correction and therefore also the SFR (in fact $E(B-V)_{\rm
cont}=0$ for both, so no extinction corrections were applied). In addition,
Q2206-BX151 is the only object that has not been observed with MOSFIRE. The
$\rm H\alpha$ flux measurement from Keck-NIRSPEC is reported by Kulas et al.
(2012), and has significant systematic uncertainties due to slit losses and
the difficulties of accurate flux calibration (Erb et al. 2006 estimated a
typical factor of $\sim 2$ slit loss correction for NIRSPEC observations of
$\rm H\alpha$ emission at $z\sim 2$).

From the stellar masses and SFRs we calculate the specific star formation
rate, sSFR $\equiv$ SFR/$M_{\star}$, finding a sample median of 10.4 Gyr-1,
more than a factor of four larger than the KBSS-MOSFIRE sample median of 2.4
Gyr-1 (Strom et al. 2017). In other words, most of the galaxies in this sample
lie significantly above the $z\sim 2$ SFR-stellar mass relation (Reddy et al.
2012; Whitaker et al. 2014, but note that samples remain incomplete at the low
masses characteristic of our targets).

The [O III]/[O II] ratios222O32 $\equiv$ [O III]$\lambda\lambda$4960,5008/[O
II]$\lambda\lambda$3727,3729, corrected for extinction. provide an estimate of
the degree of excitation in the H II regions. As expected given the sample
selection criterion of high [O III]/$\rm H\beta$ ratios, our targets fall at
the upper end of the O32 distribution for the KBSS sample (Strom et al. 2017).
High [O III]/$\rm H\beta$ and O32 are both associated with high equivalent
width $\rm Ly\alpha$ emission (e.g. Nakajima et al. 2016; Trainor et al.
2019).

We calculate the electron density from the [O II]$\lambda$3729/[O
II]$\lambda$3727 ratio, which has a median value of 1.01, slightly lower than
the ratios of 1.13–1.16 found for composite spectra of KBSS galaxies at $z\sim
2$ (Steidel et al. 2014; Strom et al. 2017). This lower ratio indicates that
many of the galaxies in our sample have higher than average electron
densities. Using the [O II] electron density calibration of Sanders et al.
(2016), we find median $n_{e}\sim 450$ cm-3. The nebular line ratios and
electron temperatures are listed in Table 2.

In summary, the galaxies studied here are relatively low-mass, highly star-
forming, and luminous, with high specific star formation rates and nebular
line ratios that place them at the upper end of their parent sample in
ionization and electron density. They are not typical of star-forming galaxies
at $z\sim 2$, but may more closely resemble galaxies observed in the
reionization era (e.g. Stark et al. 2017).

###  2.2 Observations

The 12 targets were observed with KCWI over the course of a number of
observing runs between September 2018 and August 2020. We used the Medium IFU
with the BL grating, which provides a field of view of 16$\farcs$5 $\times$
20$\farcs$4 and spectral resolution $R\approx 1800$. As detailed in Table 1,
total integration times were approximately five hours per target, divided into
individual 1200 s exposures between which we rotated the field by
$10$–$90^{\circ}$.

###  2.3 Data reduction

The KCWI data were reduced using procedures described in detail by Chen et al.
(2021), but we give an overview of the method here. Each KCWI exposure was
reduced using the official data reduction pipeline (DRP) written in
IDL.333https://github.com/Keck-DataReductionPipelines/KcwiDRP The DRP conducts
overscan and bias subtraction, cosmic-ray removal, flat-fielding, sky
subtraction, differential atmospheric refraction correction, and flux
calibration, and assembles 2D spectra of the slices into a 3D data cube. A
median-filtered cube was constructed for each data cube using a running boxcar
filter of size 0$\farcs$69 $\times$ 4$\farcs$6 $\times$ 100 Å. This median-
filtered cube was subtracted from the original cube to remove low-frequency
scattered light in both the spatial and spectral dimensions.

The world coordinate system (WCS) of the post-DRP data cubes was corrected by
cross-correlating the pseudo-white-light images of the data cubes with each
other. Data cubes of multiple exposures for the same target were rotated to
the north-up direction and resampled onto a common 3D grid of 0$\farcs$3
$\times$ 0$\farcs$3 $\times$ 1 Å. The resampling was conducted using the
“drizzle” method in the Montage package,444http://montage.ipac.caltech.edu/
with a drizzle factor of 0.7. Finally, individually resampled data cubes were
weighted by exposure time and averaged, creating the final data cube for each
target.

=30mm

Table 3: $\rm Ly\alpha$ Measurements

ID | $F^{\rm tot}_{\rm Ly\alpha}$ | Log ($L^{\rm tot}_{\rm Ly\alpha}$) |  $\rm Ly\alpha$/$\rm H\alpha$ | ($\rm Ly\alpha$/$\rm H\alpha$)int | $f_{\rm esc}^{\rm Ly\alpha}$ | $r^{\rm eff}_{\rm Ly\alpha}$ | $(F_{\rm blue}/F_{\rm red})^{\rm 1D}$ | $\Delta v^{\rm 1D}_{\rm peak}$ | $f^{\rm 1D}_{\rm tr}$ | $W^{\rm 1D}_{\rm Ly\alpha}$ | $W^{\rm tot}_{\rm Ly\alpha}$
---|---|---|---|---|---|---|---|---|---|---|---
| ($10^{-16}$ $\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}$) | (erg s-1) |  |  |  | (kpc) |  | (km s-1) |  | (Å) | (Å)
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12)
Q0142-BX165 | 3.69 $\pm$ 0.03 | 43.2 | 1.57 $\pm$ 0.22 | 8.4 | 0.19 $\pm$ 0.027 | 22.6 | 54.1 $\pm$ 1.9 | 84.5 $\pm$ 2.7 | 0.17 $\pm$ 0.005 | 658 $\pm$ 10 | 0.021 $\pm$ 0.002
Q0142-BX186 | 0.36 $\pm$ 0.02 | 42.2 | 0.31 $\pm$ 0.04 | 8.4 | 0.04 $\pm$ 0.005 | 16.4 | 30.9 $\pm$ 6.0 | 49.6 $\pm$ 8.0 | 0.32 $\pm$ 0.067 | 687 $\pm$ 31 | 0.019 $\pm$ 0.016
Q0207-BX87 | 3.72 $\pm$ 0.03 | 43.1 | 3.88 $\pm$ 0.69 | 9.0 | 0.43 $\pm$ 0.076 | 25.6 | 73.3 $\pm$ 4.1 | 130.3 $\pm$ 6.9 | 0.23 $\pm$ 0.007 | 522 $\pm$ 13 | 0.074 $\pm$ 0.002
Q0207-BX144 | 4.61 $\pm$ 0.04 | 43.2 | 2.06 $\pm$ 0.22 | 8.5 | 0.24 $\pm$ 0.026 | 27.1 | 44.1 $\pm$ 2.5 | 110.2 $\pm$ 5.1 | 0.18 $\pm$ 0.012 | 631 $\pm$ 7 |  $-$0.002 $\pm$ 0.003
Q0449-BX110 | 1.98 $\pm$ 0.03 | 42.9 | 1.10 $\pm$ 0.09 | 9.0 | 0.12 $\pm$ 0.010 | 20.9 | 49.8 $\pm$ 4.0 | 87.6 $\pm$ 5.4 | 0.42 $\pm$ 0.021 | 700 $\pm$ 13 | 0.038 $\pm$ 0.004
Q0449-BX115 | 1.73 $\pm$ 0.03 | 42.9 | 8.57 $\pm$ 2.86aaCalculated from $\rm H\beta$ flux, assuming $E(B-V)=0.0$ from SED fitting and an intrinsic ratio $\rm H\alpha/\rm H\beta=2.79$.  | 8.9 | 0.96 $\pm$ 0.320 | 17.1 | 98.9 $\pm$ 10.6 | 129.7 $\pm$ 12.1 | 0.21 $\pm$ 0.011 | 518 $\pm$ 8 | 0.029 $\pm$ 0.003
Q0821-MD36 | 2.11 $\pm$ 0.03 | 43.1 | 1.13 $\pm$ 0.20 | 9.1 | 0.12 $\pm$ 0.022 | 20.2 | 86.8 $\pm$ 7.9 | 145.8 $\pm$ 12.3 | 0.15 $\pm$ 0.009 | 352 $\pm$ 5 | 0.091 $\pm$ 0.003
Q1549-BX102 | 1.87 $\pm$ 0.03 | 42.9 | 2.71 $\pm$ 1.13 | 8.9 | 0.30 $\pm$ 0.127 | 21.6 | 44.9 $\pm$ 3.8 | 78.0 $\pm$ 5.2 | 0.43 $\pm$ 0.019 | 649 $\pm$ 9 | 0.016 $\pm$ 0.004
Q1700-BX729 | 1.19 $\pm$ 0.02 | 42.7 | 0.53 $\pm$ 0.06 | 8.9 | 0.06 $\pm$ 0.007 | 20.6 | 24.6 $\pm$ 2.0 | 43.6 $\pm$ 2.9 | 0.11 $\pm$ 0.017 | 853 $\pm$ 29 |  $-$0.022 $\pm$ 0.007
Q2206-BX151 | 6.92 $\pm$ 0.05 | 43.4 | 10.98 $\pm$ 0.86 | 8.9 | 1.23 $\pm$ 0.096 | 29.8 | 122.5 $\pm$ 8.1 | 195.0 $\pm$ 13.1 | 0.31 $\pm$ 0.007 | 598 $\pm$ 4 | 0.029 $\pm$ 0.002
Q2343-BX418 | 3.97 $\pm$ 0.03 | 43.2 | 2.71 $\pm$ 0.24 | 9.1 | 0.30 $\pm$ 0.026 | 19.8 | 64.3 $\pm$ 2.2 | 80.6 $\pm$ 2.3 | 0.38 $\pm$ 0.009 | 563 $\pm$ 2 | 0.029 $\pm$ 0.001
Q2343-BX660 | 2.51 $\pm$ 0.03 | 43.0 | 1.54 $\pm$ 0.17 | 8.7 | 0.18 $\pm$ 0.020 | 26.0 | 33.5 $\pm$ 1.7 | 110.4 $\pm$ 4.7 | 0.09 $\pm$ 0.011 | 690 $\pm$ 11 |  $-$0.016 $\pm$ 0.003

††footnotetext: Columns: (1) Galaxy ID; (2) Total $\rm Ly\alpha$ flux in halo;
(3) Total $\rm Ly\alpha$ luminosity of halo; (4) Observed $\rm Ly\alpha$/$\rm
H\alpha$ ratio, with $\rm H\alpha$ corrected for extinction; (5) Intrinsic
$\rm Ly\alpha$/$\rm H\alpha$ ratio calculated from electron density; (6) $\rm
Ly\alpha$ escape fraction; (7) Effective circular radius of halo, calculated
from total area $A$ as $r^{\rm eff}_{\rm Ly\alpha}=(A/\pi)^{1/2}$; (8)
Blue/red peak flux ratio of optimally-extracted 1D spectrum; (9) Peak
separation of optimally-extracted 1D spectrum; (10) Fraction of the total $\rm
Ly\alpha$ emission within $\pm 100$ km s-1 of the trough between the peaks;
(11) Rest-frame $\rm Ly\alpha$ equivalent width from optimally-extracted 1D
spectrum; (12) Total rest-frame equivalent width of $\rm Ly\alpha$ halo.

##  3 Global $\rm Ly\alpha$ measurements

In this section we describe the global properties of the $\rm Ly\alpha$
emission measured by KCWI. We begin with one-dimensional spectra designed to
optimize the continuum S/N, for comparison with single slit studies. We define
isophotal apertures by running SExtractor (Bertin & Arnouts 1996) in detection
mode on the collapsed, white-light images, and then extract spectra from these
apertures, weighting by (S/N)2. The resulting spectra are then rescaled to
match the total aperture flux. Circularized radii of the apertures range from
0$\farcs$9 to 2$\farcs$2, with all but three within 0$\farcs$3 of the sample
median of 1$\farcs$6\. At the median redshift of the sample, these spectra
cover rest-frame wavelengths $\sim$1060–1660 Å, and include a number of
interstellar absorption lines which we model along with the $\rm Ly\alpha$
emission in Section 6.1.

We show the continuum-subtracted $\rm Ly\alpha$ profiles from these spectra in
Figure 1, which demonstrates that all objects in the sample have double-peaked
profiles with a dominant red peak. The $\rm Ly\alpha$-adjacent continuum is
defined as the median flux density in two windows on either side of the line,
spanning 1199–1210 Å ($-4000$ to $-1400$ km s-1) on the blue side and
1225–1236 Å ($+2300$ to $+5000$ km s-1) on the red side. Given the generally
high continuum S/N of the optimally extracted spectra and the lack of
underlying absorption, these relatively narrow windows provide an effective
measurement of the continuum around the line, as can be seen by assessing the
continuum subtraction in Figure 1.

In order to measure equivalent widths, we integrate the line between the
limits at which it reaches the continuum and divide the resulting flux by the
continuum level determined above. The resulting rest-frame $\rm Ly\alpha$
equivalent widths are given in the column labeled $W^{\rm 1D}_{\rm Ly\alpha}$
in Table 3. As previously known and by design, all are above the canonical
$\rm Ly\alpha$-emitter threshold of $W_{\rm Ly\alpha}>20$ Å.

We next create pseudo-narrowband continuum-subtracted $\rm Ly\alpha$ surface
brightness images for each object in the sample. We first identify the spatial
peak of the $\rm Ly\alpha$ emission in each data cube, and then extract the
summed one-dimensional $\rm Ly\alpha$ profile of a large (2$\farcs$4 in
diameter) region centered on this peak. We measure the wavelengths at which
the $\rm Ly\alpha$ emission from this large region meets the continuum on
either side of the line, typically $\sim-900$ to $+1200$ km s-1, and use these
as the wavelength limits of a $10\farcs 5\times 10\farcs 5$ (i.e. $35\times
35$ 0$\farcs$3 pixels) subcube centered on the $\rm Ly\alpha$ peak. We also
extract blue and red subcubes with the same spatial size and spectral widths
of 20 Å in the rest frame from either side of the $\rm Ly\alpha$ emission,
from which we measure the continuum level (because we are here measuring the
continuum of individual spaxels rather than a spatially integrated region as
we did above, we use slightly wider windows to increase the S/N of the
measurement). The median of the blue and red subcubes along the wavelength
axis results in a continuum image, which we subtract from the $\rm Ly\alpha$
subcube to create an emission-only cube. Finally, this cube is integrated
along the wavelength axis to construct the continuum-subtracted $\rm Ly\alpha$
surface brightness images shown in Figure 2. We note that, in general,
accurate modeling and subtraction of the continuum underlying $\rm Ly\alpha$
emission can be challenging due to the complex nature of the line profiles,
which often display a superposition of emission and absorption (e.g. Shapley
et al. 2003; Kornei et al. 2010). However, the targets in the current sample
have simpler profiles with strong emission and no detectable absorption,
enabling effective continuum subtraction with the simple method described
here.

The $\rm Ly\alpha$ emission shown in Figure 2 is significantly more extended
than the underlying UV continuum in all cases, as can be seen by comparing the
white ($\rm Ly\alpha$) and black (continuum) contours in Figure 2. This
comparison also shows that in most cases the $\rm Ly\alpha$ and continuum
peaks are spatially coincident. We measure the total $\rm Ly\alpha$ flux of
each extended halo by summing the largest connected region with S/N $>1$
(corresponding to a surface brightness of $\sim 1$–$2\times 10^{-18}$
$\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{arcsec}^{-2}$ for
targets with approximately 5 hrs of integration), and calculate an effective
circular radius $r^{\rm eff}_{\rm Ly\alpha}=(A/\pi)^{1/2}$, where $A$ is the
area of the region. The total $\rm Ly\alpha$ fluxes range from 0.4 $7\times
10^{-16}$ $\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}$, corresponding to
$\rm Ly\alpha$ luminosities ranging from $1.6\times 10^{42}$ to $2.6\times
10^{43}$ erg s-1 with median $10^{43}$ erg s-1, while the radii vary between
16 and 30 kpc.

$\rm Ly\alpha$ escape fractions are computed by comparing the total $\rm
Ly\alpha$ flux with the predicted flux calculated from the dust-corrected $\rm
H\alpha$ emission555We use $\rm H\beta$ emission for Q0449-BX115, assuming
$\rm H\alpha$/$\rm H\beta$ =2.79. and the theoretical $\rm Ly\alpha$/$\rm
H\alpha$ ratio. Estimates of the $\rm Ly\alpha$ escape fraction typically
assume an intrinsic ratio $(\rm Ly\alpha/\rm H\alpha)_{\rm int}=8.7$ (e.g.
Hayes 2015), but the $\rm Ly\alpha$/$\rm H\alpha$ ratio is density-dependent,
increasing at higher electron densities as collisional processes suppress two-
photon continuum emission. For the range of densities in our sample,
$n_{e}\sim 100$ – 700 cm-3, $(\rm Ly\alpha/\rm H\alpha)_{\rm int}$ ranges from
8.4 to 9.1.666We assume $T_{e}=15,000$ K. The dependence of the $\rm
Ly\alpha$/$\rm H\alpha$ ratio on temperature is smaller than the dependence on
density, $\sim 2$% between 10,000 and 20,000 K. We therefore adopt an
intrinsic ratio for each object that depends on the electron density; these
ratios are listed in the column labeled $(\rm Ly\alpha/\rm H\alpha)_{\rm int}$
in Table 3. For the three galaxies that do not have measurements of $n_{e}$
(Q0449-BX115, Q1700-BX729 and Q2206-BX151) we adopt the sample median of
$n_{e}\sim 450$ cm-3, corresponding to $(\rm Ly\alpha/\rm H\alpha)_{\rm
int}=8.9$.

The resulting escape fractions range from 0.04 to 1.2, with a median of 0.22.
We measure $f_{\rm esc}^{\rm Ly\alpha}=1.23\pm 0.096$ for Q2206-BX151,
suggesting $\sim 100$% of the $\rm Ly\alpha$ emission escapes; however, as
discussed in Section 2.1, the $\rm H\alpha$ flux for this object has
significant systematic uncertainties and is likely underestimated, leading to
an overestimate of the escape fraction. Nevertheless, Q2206-BX151 does have
the largest $\rm Ly\alpha$ luminosity, halo size, and equivalent width in the
sample, so a high escape fraction may also be expected. Similarly, the second
highest escape fraction measured, $f_{\rm esc}^{\rm Ly\alpha}=0.96\pm 0.32$
for Q0449-BX115, may also be overestimated due to the use of $E(B-V)_{\rm
cont}$ rather than $E(B-V)_{\rm neb}$ for the dust correction; in practice,
$E(B-V)_{\rm cont}=0$ for Q0449-BX115, so no dust correction was applied.

We next calculate the total $\rm Ly\alpha$ equivalent width by dividing the
total $\rm Ly\alpha$ flux by the median continuum flux density of the
optimally extracted one-dimensional spectra, determined as described above.
After conversion to the rest frame, these total equivalent widths are larger
than the 1D equivalent widths measured above by factors ranging from 1.3 to
3.3, with a median of 1.7. The largest total equivalent width in the sample is
$\sim 200$ Å (for Q2206-BX151), roughly the maximum value expected from a
normal stellar population (Charlot & Fall 1993). All $\rm Ly\alpha$
measurements are given in Table 3.

Finally, we note that in this work we focus on the spectral properties of the
extended $\rm Ly\alpha$ emission in this strongly $\rm Ly\alpha$-emitting
subset of the KBSS sample; an analysis of the structural and spectral
properties of $\rm Ly\alpha$ emission in the full KBSS-KCWI sample will be
presented elsewhere.

Figure 3: Representative results of the Voronoi binning of the $\rm Ly\alpha$
surface brightness images. The bins include only individual pixels with S/N
$\geq 1$, and are constructed to achieve a target S/N $\geq 3$. Bins generally
consist of single pixels, except in the outer halo where each object has a few
bins of $\sim 2$–5 pixels. Contours showing $\rm Ly\alpha$ surface brightness
are the same as in Figure 2, and show that regions reaching the target S/N
generally have surface brightness $\gtrsim 5\times 10^{-18}$
$\mathrm{erg}\,\mathrm{s}^{-1}\,\mathrm{cm}^{-2}\,\mathrm{arcsec}^{-2}$
(dashed contour).

##  4 The spectral profiles of spatially extended $\rm Ly\alpha$ emission

In this section we analyze the spatial variations of the spectral $\rm
Ly\alpha$ profiles across the extended halos using individual spaxels and
small bins, focusing on the flux ratio and velocity separation of the two
peaks. We first apply an adaptive two-dimensional Voronoi binning procedure to
increase the S/N in the outer regions of the halos using the python package
vorbin, an implementation of the method described in detail by Cappellari &
Copin (2003). Beginning with the highest S/N pixel, this routine works its way
outward to lower S/N regions until a pixel with S/N lower than a specified
target threshold is reached; a bin is then constructed from adjacent pixels,
seeking to match the target S/N. We run the binning on the $\rm Ly\alpha$
surface brightness images described in Section 3 above, avoiding the inclusion
of pure noise by using only individual pixels with S/N $>1$ and setting a
target for each bin of S/N $\geq 3$. In practice, most ($\sim 90$%) of the
resulting bins consist of single pixels, and very few (0–3 per object) contain
four or more pixels. Three representative examples of the results of the
Voronoi binning are shown in Figure 3, with only bins with S/N $\geq 3$ shown
and each bin indicated by a different color.

Figure 4: The $\rm Ly\alpha$ profiles (black) and asymmetric Gaussian fits
(orange dashed lines) for individual spaxels and Voronoi bins in the NW
portion of the halo of Q0449-BX110. The two spaxels in the grey box in the top
row comprise a single bin, and the two spectra plotted within the box are
identical. The blue shaded regions indicate that the blue side of the line is
undetected. The image at upper right shows the $\rm Ly\alpha$ surface
brightness of the full modeled region of Q0449-BX110 on a logarithmic scale,
with the portion shown here outlined in orange. The $\rm Ly\alpha$ profiles
are shown in the order of their corresponding spaxels in the orange box.
Contours of $\rm Ly\alpha$ surface brightness are the same as in Figure 2.

Once the bins are identified, we extract the spectrum and variance of each bin
from the continuum-subtracted $\rm Ly\alpha$ datacube, using the mean of the
individual spaxels for bins consisting of more than one spaxel. We then fit
the $\rm Ly\alpha$ profile of each bin with a double asymmetric Gaussian
function in velocity space, defined as

$\displaystyle f(v)={}A_{\rm blue}$ $\displaystyle\exp\left(\frac{-(v-v_{\rm
0,blue})^{2}}{2\sigma_{\rm blue}^{2}}\right)+$ (1) $\displaystyle A_{\rm red}$
$\displaystyle\exp\left(\frac{-(v-v_{\rm 0,red})^{2}}{2\sigma_{\rm
red}^{2}}\right),$

where $A_{\rm blue}$ and $A_{\rm red}$ and $v_{\rm 0,blue}$ and $v_{\rm
0,red}$ are the amplitudes and peak velocities of the blue and red components
respectively. The asymmetric line width $\sigma$ is further defined as
$\sigma=a(v-v_{0})+d$, where $a$ and $d$ describe the asymmetry and width of
the profile respectively. A single asymmetric Gaussian has been previously
used to fit $\rm Ly\alpha$ profiles at high redshift by Shibuya et al. (2014)
and Leclercq et al. (2020); here we introduce separate components for the blue
and red parts of the lines, given the strong double-peaked nature of our
sources. Double asymmetric Gaussian fits to the spatially integrated $\rm
Ly\alpha$ spectra are shown in Figure 1, demonstrating that they generally
provide an excellent representation of the line profile.

For each spaxel we measure the flux ratio of the blue and red peaks by
integrating each side of the line between the trough between the peaks (or
zero velocity, in the rare cases in which the trough is not present) and the
point at which the S/N drops below unity. For spaxels with significant
detections of both peaks, we use the parameters of the fit to define the $\rm
Ly\alpha$ peak separation for each spaxel as $\Delta v_{\rm peak}=v_{\rm
0,red}-v_{\rm 0,blue}$. Flux uncertainties are determined via error
propagation of the variance cube, while uncertainties on the peak separation
result from the covariance matrix of the Gaussian fit. We also tested a Markov
Chain Monte Carlo approach to the fitting and uncertainties, but found that it
did not change the results significantly while being very time-consuming given
the number of fits involved.

Figure 4 shows an example of the spectra of individual spaxels and Voronoi
bins for a portion of the halo of Q0449-BX110, chosen because it is near the
sample median in both total $\rm Ly\alpha$ flux and halo size. The asymmetric
Gaussian fits to each spaxel are overplotted in orange.

Maps of the $\rm Ly\alpha$ blue/red flux ratio and peak separation resulting
from these measurements for the 12 galaxies are shown in Figures 5 and 6,
where each lettered panel shows the peak ratio at the top and the peak
separation at bottom. Higher flux ratios (i.e. those with a stronger blue
peak) and smaller peak separations are indicated in blue. Most objects have
peak ratios ranging from $\sim 0.1$–$1$ and peak separations ranging from
$\sim 300$ to $\sim 700$ km s-1, although a few have lower peak ratios
(Q0449-BX115, Q0821-MD36, Q1700-BX729) or larger separations (Q1700-BX729).
The maps also indicate that regions of higher peak ratio and lower peak
separations are generally found in the outer parts of the halos, although not
always in the same regions.

Figure 5: Maps of the $\rm Ly\alpha$ peak ratio and separation for the first
six objects in the sample. Each lettered panel (a) through (f) shows the peak
ratio at the top and the peak separation at bottom for a given object. The
blue peak is undetected for spaxels marked with an $\times$, the red peak is
undetected for spaxels marked with a $+$, and the color of these spaxels
indicates the 1$\sigma$ upper limit in the case of blue non-detections and the
1$\sigma$ lower limit in the case of red non-detections. We measure the peak
separation only for spaxels for which both peaks are detected. White and black
contours show the $\rm Ly\alpha$ and continuum surface brightnesses
respectively, as in Figure 2.

Figure 6: Maps of the $\rm Ly\alpha$ peak ratio and separation for objects
7–12. Each lettered panel (g) through (l) shows the peak ratio at the top and
the peak separation at bottom for a given object. The blue peak is undetected
for spaxels marked with an $\times$, the red peak is undetected for spaxels
marked with a $+$, and the color of these spaxels indicates the 1$\sigma$
upper limit in the case of blue non-detections and the 1$\sigma$ lower limit
in the case of red non-detections. We measure the peak separation only for
spaxels for which both peaks are detected. White and black contours show the
$\rm Ly\alpha$ and continuum surface brightnesses respectively, as in Figure
2.

Figure 7: $\rm Ly\alpha$ peak ratio vs. separation. Each point represents a
single spaxel or Voronoi bin, and is color-coded by its distance from the peak
of the $\rm Ly\alpha$ surface brightness image. Spearman $r$ and $p$ values
are given in the upper right of each panel, with significant ($>3\sigma$)
correlations labeled in red. With $p=2.83\times 10^{-3}$, Q1700-BX729 has
2.99$\sigma$ significance.

Figure 8: Significance levels from pairwise Spearman correlation tests of
distance from the peak of the $\rm Ly\alpha$ surface brightness image, $\rm
Ly\alpha$ surface brightness, peak ratio and peak separation.

We study the relationship between the peak ratio and separation further in
Figure 7, in which we plot the ratio vs. separation for each spaxel or bin of
each of the 12 objects, with each point color-coded by its distance from the
center, defined as the $\rm Ly\alpha$ surface brightness peak. This color-
coding again shows that higher ratios and lower separations tend to be found
at larger radii. For each object we measure the Spearman correlation
coefficient $r$ and the probability $p$ of the null hypothesis that the peak
ratio and separation are uncorrelated; these results are reported in the upper
right of each panel, with significant correlations ($>3\sigma$, $p<0.0027$)
labeled in red. Five of the twelve objects in the sample show significant
anti-correlations, such that a higher peak ratio is associated with a lower
peak separation; one additional object, Q1700-BX729, shows a 2.99$\sigma$
correlation.

To gain further insight into the connection between peak ratio and separation,
we perform pairwise Spearman correlation tests on the four quantities we
measure for each spaxel or Voronoi bin: distance from the center, $\rm
Ly\alpha$ surface brightness, peak ratio, and peak separation. This results in
six correlation coefficients for each object, and we show the significance
levels of these correlations in Figure 8. In addition to the expected
extremely strong correlation between distance and surface brightness, we see
strong ($>5\sigma$) correlations between peak ratio and distance or surface
brightness for most (8/12) of the sample. Correlations between peak separation
and distance or surface brightness are usually weaker, but are present at
$>3\sigma$ in 8/12 objects. For two galaxies in the sample (Q0449-BX110 and
Q1549-BX102) the separation is more strongly correlated with distance and
surface brightness than the ratio. All objects except Q0142-BX186 have at
least one $>3\sigma$ correlation in addition to that of distance and surface
brightness, and there is no obvious preference for either distance or surface
brightness to be more strongly correlated with the line profile properties.
The galaxy showing no correlations with $\rm Ly\alpha$ peak properties,
Q0142-BX186, is the faintest object in the sample, with total $\rm Ly\alpha$
flux more than three times lower than the second faintest (Q1700-BX729),
demonstrating that high S/N over a relatively large area is needed to detect
these trends.

These results suggest that the correlations between peak ratio and separation
are largely driven by the underlying tendencies of these quantities to
increase and decrease respectively in the outer, fainter parts of the halos.
We will discuss the relationship between the ratio and separation in more
detail in the following sections.

Figure 9: Normalized annular $\rm Ly\alpha$ profiles, constructed by binning
all spaxels that have S/N $>2$ in the continuum-subtracted $\rm Ly\alpha$
images in single spaxel (0$\farcs$3) radial increments. The spectra are color-
coded by radius, with the inner portions of the halo in red and the outer
portions in blue. The legend in each panel gives the median radius of each
bin. For most of the sample, the blue-to-red peak ratio increases and the
depth of the trough between the peaks decreases with increasing radius.

##  5 Spatially averaged $\rm Ly\alpha$ profiles

In order to determine general trends and study the spectral properties of the
$\rm Ly\alpha$ halos to larger radii than can be measured with individual
spaxels and the small Voronoi bins, we also construct binned spectra of larger
regions, using both the entire halo and smaller regions chosen based on their
spectral properties.

###  5.1 Annular Ly$\alpha$ profiles

We first study the average variation of the $\rm Ly\alpha$ profile as a
function of radius by making annular spectra binned by radius for all objects
in the sample. Beginning with the central, highest surface brightness spaxel
and including all spaxels with S/N $>2$ in the continuum-subtracted $\rm
Ly\alpha$ images, we bin each halo in single spaxel (0$\farcs$3) radial
increments. The spectra of all the spaxels in each bin are then summed, and
the resulting $\rm Ly\alpha$ profile is normalized to a total flux of 1. This
normalization enables a straightforward visual examination of changes in the
shape of the profile with radius, and is also used to format the spectra for
the radiative transfer modeling discussed in Section 6. The normalized spectra
are shown in Figure 9, color-coded by radius with the central portions of each
halo in red and the outer portions in blue.

The increasing strength of the blue peak relative to the red peak with
increasing radius is clearly apparent for most of the objects in the sample.
It is also clear that the depth of the trough between the two peaks decreases
with increasing radius for most of the sample. Although generally less obvious
to the eye, the trend of decreasing peak separation with increasing radius is
also apparent in many of the sources. We quantify these trends by fitting
double asymmetric Gaussian profiles to the binned annular spectra as described
in Section 4 above, measuring the average peak ratio and separation as a
function of both radius and average surface brightness. We also quantify the
depth of the trough between the peaks by measuring $f_{\rm tr}$, defined as
the fraction of the total emission within $\pm 100$ km s-1 of the
trough.777$f_{\rm tr}$ differs slightly from the quantity $f_{\rm cen}$
defined by Naidu et al. (2022), who measure the fraction of flux escaping
within $\pm 100$ km s-1 of the systemic velocity; we instead measure the flux
on either side of the trough to account for the fact that the trough is
occasionally slightly offset from zero velocity (e.g. Q2343-BX660).

The results are shown for all objects in Figure 10, and generally confirm
expectations from visual inspection of the spectra. Central blue-to-red flux
ratios are $\sim 0.2\pm 0.2$, and the average ratio increases consistently
with radius for most of the galaxies in the sample; all objects that can be
measured at a radius beyond $\sim 16$ kpc have flux ratios $>0.6$ at that
radius. The trough flux fraction $f_{\rm tr}$ ranges from $<0$ to $\sim 0.1$
at the center of the halos, and rises consistently with radius for most
objects. The largest measured values are $f_{\rm tr}\sim 0.2$, found in the
outer halos of Q0207-BX87 and Q2343-BX418.

Trends with the average peak separation are somewhat more complicated. Most
(10/12) of the halos have a central peak separation of $\sim 500$–700 km s-1,
with the exceptions of Q0821-MD36 (365 km s-1) and Q1700-BX729 (835 km s-1).
In most cases the average peak separation decreases with radius, with a
typical change of $\sim-100$ km s-1 such that peak separations in the outer
halo are $\sim 400$–600 km s-1; however, a few objects (e.g. Q0142-BX165 and
Q2343-BX418) show steeper gradients. Two galaxies in the sample (Q2206-BX151
and Q2343-BX660) also show an increase in the peak separation at the largest
radius; in both cases these increases are due to small regions with large
separations at large radius, as can be seen in Figure 6. Unsuprisingly,
however, the peak separation is closely related to the trough depth $f_{\rm
tr}$, decreasing as $f_{\rm tr}$ increases.

Figure 10: Results of the line profile measurements of the annular spectra
described in Section 5.1. Top row: $\rm Ly\alpha$ peak ratio vs. radius; $\rm
Ly\alpha$ peak separation vs. radius; and peak separation vs. ratio. Middle
row: $\rm Ly\alpha$ peak ratio vs. normalized $\rm Ly\alpha$ surface
brightness; $\rm Ly\alpha$ peak separation vs. surface brightness; and surface
brightness vs. radius. Bottom row: The fraction of total flux within $\pm 100$
km s-1 of the trough between the peaks $f_{\rm tr}$ vs. radius; $f_{\rm tr}$
vs. peak separation.

Figure 11: Top two rows: Blue-to-red $\rm Ly\alpha$ flux ratios from radially
binned spectra of the 60∘ angular regions that maximize (dark blue circles)
and minimize (light blue triangles) the gradient in peak ratio with radius.
Second two rows: Same as first two rows, for measurements of the $\rm
Ly\alpha$ peak separation, with maximum gradients indicated by gold squares
and minimum gradients by light yellow diamonds. The annular averages from
Figure 10 are also shown as red stars in all panels. See Section 5.2 for
details.

###  5.2 Gradients in Ly$\alpha$ peak ratio and separation

While the binned, annular profiles described above are useful to characterize
general trends in the extended $\rm Ly\alpha$ emission, they also wash out the
spectral variations seen in different parts of individual halos. As is readily
apparent from the maps of peak ratio and separation in Figures 5 and 6, the
$\rm Ly\alpha$ profiles across the halos are not radially symmetric, and there
are significant differences in both peak ratio and separation at different
position angles in a given halo. We therefore characterize the variations in
the $\rm Ly\alpha$ profile within individual halos by binning smaller regions,
using seven of the eight brightest sources in the sample (we do not include
Q0821-MD36, for which the blue peak is too weak to obtain useful measurements
from binning smaller regions).

Our goal is to construct a series of binned spectra that maximize or minimize
the gradients in peak ratio or separation from the center to the outskirts of
the halo. Again beginning with all spaxels with S/N $>2$ in the $\rm Ly\alpha$
images, we then take a subset of each halo corresponding to a 60∘ angular
region (chosen to encompass a large enough region to increase the S/N by
binning while still isolating different parts of the halos). As with the
annular spectra, we radially bin the datacube in this region in single spaxel
annular increments and measure the peak ratio and separation of each of the
resulting $\rm Ly\alpha$ profiles. We then rotate the 60∘ region by 10∘ and
repeat the process until the entire halo has been covered.

We next measure the peak ratio and separation for each of the resulting 36
spectra, and locate the regions of maximum and minimum gradients in peak ratio
and separation by identifying the two regions for which the difference in peak
ratio with radius is maximized, and the two regions for which the difference
in separation is maximized. For the peak ratio, the maximum gradient
corresponds to the largest increase from the center to the outskirts, while
for the separation it is the largest decrease. In other words, the steepest
peak ratio gradient is found in the direction of the highest blue-to-red flux
ratio, and the steepest peak separation gradient is found in the direction of
the narrowest peak separation.

The results of this process are shown in Figure 11, in which we plot the
maximum and minimum ratio gradients in the top two rows and the maximum and
minimum separation gradients in the bottom two rows, along with the annular
averages from Figure 10. Although the sample for which these measurements are
feasible is small, this exercise shows that all of the halos have a region for
which the peak ratio increases with radius, and a region for which the
separation decreases with radius. Notably, however, in all cases the angular
regions corresponding to these two maximum gradients do not overlap; this
result is consistent with the finding in Section 4 above that the correlation
between peak ratio and separation is largely due to the underlying
relationship of both with radius.

Turning to the minimum gradients, most of the halos also have at least one
sightline for which the increase in peak ratio with radius is small or
nonexistent, and at least one sightline for which the peak separation is
relatively flat with radius (or even rising, in the case of Q0207-BX144).
Unlike the maximum gradients, there is some overlap between the regions of
minimum gradient; for four of the seven objects, the minumum gradient regions
overlap by 10–30∘. The minimum gradients show that most halos have regions for
which the $\rm Ly\alpha$ line profile does not follow the average trends. We
discuss the implications of this observation further in Section 7, informed by
the results of spatially resolved modeling of the $\rm Ly\alpha$ emission.

##  6 Modeling Ly$\alpha$ emission and low ionization interstellar absorption
lines

In the previous sections we have shown that the spectral morphology of $\rm
Ly\alpha$ emission changes significantly across the extended halos. On
average, the blue-to-red peak flux ratio increases, the peak separation
decreases, and the fraction of the total flux emerging between the two peaks
increases with increasing radius; there are, however, variations in these
patterns with azimuthal angle within a given halo. In this section we further
examine both the spatially resolved $\rm Ly\alpha$ profiles and the “down-the-
barrel” rest-frame UV low ionization interstellar metal absorption lines using
physical models. This analysis will help us construct a consistent picture of
the ISM and CGM of the galaxies in our sample.

###  6.1 Ly$\alpha$ radiative transfer modeling

To extract physical properties of the gas in the halos from the observed $\rm
Ly\alpha$ profiles, we perform Monte Carlo radiative transfer (MCRT) modeling
of the $\rm Ly\alpha$ line. In contrast to the majority of previous studies in
which spatially integrated $\rm Ly\alpha$ spectra are modeled, in this work we
attempt to fully leverage the power of KCWI and reproduce the spatially
varying trends of the observed $\rm Ly\alpha$ profiles.

Following a similar methodology to Li et al. (2022), we model the spatially
resolved $\rm Ly\alpha$ profiles using the multiphase, clumpy model. Each
model is a 3D spherically symmetric region that emulates a galactic halo with
a $\rm Ly\alpha$ emitting source located at its center and two phases of gas:
cool ($\sim 10^{4}$ K) H I clumps and a hot ($\sim 10^{6}$ K), highly ionized
inter-clump medium (ICM). As we will show below, such a hot, diffuse, low-
density H I component is necessary to reproduce the observed $\rm Ly\alpha$
profiles, primarily by producing additional absorption near the $\rm Ly\alpha$
line center.888In the multiphase, clumpy model, the flux at line center of the
emergent $\rm Ly\alpha$ spectra is predominantly controlled by the residual
${\rm H\,{\textsc{i}}}$ number density in the static ICM component. Without
the ICM, a significant number of photons will escape at the line center. In
reality, such a component may correspond to the low column density absorbers
(${\rm log}\,N_{\rm HI}\lesssim 10^{17}\rm cm^{-2}$) that provide additional
$\rm Ly\alpha$ scatterings in a galactic outflow (see e.g. Section 7.3 of
Dijkstra & Kramer 2012). After interacting with these two phases of gas, $\rm
Ly\alpha$ photons that escape from different impact parameters can be
separated into different spatial bins and the emergent spectra can then be
compared to the corresponding observed spatially resolved $\rm Ly\alpha$
profiles.

In practice, we construct a grid of multiphase, clumpy models for fitting the
$\rm Ly\alpha$ spectra by varying the five most important physical parameters:
$F_{\rm V}$, the volume filling factor of the clumps; ${\rm log}\,N_{\rm
HI,\,{\rm cl}}$, the H I column density of the clumps; $\sigma_{\rm cl}$, the
random velocity dispersion of the clumps; $v_{\rm cl}$, the radial outflow
velocity of the clumps; and ${\rm log}\,n_{\rm HI,\,{\rm ICM}},$ the residual
H I number density in the ICM.999We have also experimented with a varying
$v_{\rm ICM}$ (the radial outflow velocity of the ICM), but found that in
almost all cases $v_{\rm ICM}\simeq 0$ is preferred by the fitting. This is
due to the prominent trough in most of the $\rm Ly\alpha$ profiles that
requires significant absorption at the line center. Therefore, we have fixed
$v_{\rm ICM}$ to zero to reduce the dimensionality of our model grid.  An
additional parameter, $\Delta v$, is used in post-processing to determine the
deviation between the best-fit systemic redshift of the $\rm Ly\alpha$
emitting source and the observed systemic redshift inferred from non-resonant
nebular emission lines. The parameter values of the model grid are summarized
in Table 4. Note that the range of $F_{\rm V}$ corresponds to a cloud covering
factor $f_{\rm c}^{\rm Ly\alpha}$ (the average number of clumps per line of
sight) of $\sim 1-10$, which is similar to or moderately larger than the
critical threshold $f_{\rm c,\,crit}^{\rm Ly\alpha}$. Here $f_{\rm
c,\,crit}^{\rm Ly\alpha}$ denotes the critical average number of clumps per
line of sight, above which the clumpy medium starts to transition to a
homogeneous medium. In other words, we are exploring a unique physical regime
($f_{\rm c}^{\rm Ly\alpha}\simeq f_{\rm c,\,crit}^{\rm Ly\alpha}$) where the
$\rm Ly\alpha$ RT in a multiphase, clumpy medium does not fully converge to
the homogeneous shell model (Gronke et al. 2016, 2017; Li & Gronke 2022).

Table 4: Parameter values of the multiphase, clumpy model grid.

Parameter | Definition | Values
---|---|---
(1) | (2) | (3)
$F_{\rm V}$ | Clump volume filling factor | (0.01, 0.04, 0.08, 0.16)
${\rm log}\,N_{\rm HI,\,cl}$ |  Clump H i column density |  (17, 17.5, 18, 18.5, 19) log cm-2
$\sigma_{\rm cl}$ | Clump velocity dispersion |  (0, 25, 50, …, 150) km s-1
$v_{\rm cl,\,\infty}$ | Clump asymptotic outflow velocity |  (500, 600, 700, 800, 900) km s-1
${\rm log}\,n_{\rm HI,\,{\rm ICM}}$ | ICM H I number density |  (-8, -7.5, -7, -6.5) log cm-3
$\Delta v$ |  Velocity shift relative to systemic $z$ |  [-120, 120] km s-1

††footnotetext: Notes. The parameter values of the model grid that we used for
fitting the $\rm Ly\alpha$ profiles. The columns are: (1) parameter name; (2)
parameter definition; (3) parameter values on the grid.

Previous work (Li et al. 2022) assumed constant radial outflow velocities, but
here we adopt a more physically realistic radially varying clump outflow
velocity profile. Our choice is inspired by Dijkstra & Kramer (2012), who find
that a radially varying velocity profile is able to better reproduce the
surface brightness (SB) profiles of $\rm Ly\alpha$ halos. Specifically, the
momentum equation of an H I clump can be written as (Murray et al. 2005;
Martin 2005):

$\frac{\mathrm{d}v(r)}{\mathrm{d}t}=-\frac{GM(r)}{r^{2}}+Ar^{-\alpha}$ (2)

where $r$ is the clump’s radial position, $v(r)$ is the clump radial outflow
velocity at $r$, $M(r)$ is the total gravitational mass within $r$, and $A$ is
a constant that characterizes the amplitude of the power-law acceleration
$r^{-\alpha}$. The acceleration of the clump is determined by two competing
terms on the right hand side, the first of which is due to gravitational
deceleration and the second of which is an empirical power-law acceleration
term (Steidel et al. 2010). Major acceleration mechanisms for the cool clumps
may include radiation pressure, ram pressure from a hot wind, and shock-
accelerated cosmic rays, which all may correspond to an $r^{-2}$ force (see
Chisholm et al. 2016, and note that the radiation pressure should be in the
optically thin regime). However, in reality, the clumps may suffer from extra
deceleration (and acceleration, see Gronke & Oh 2020) due to their interaction
with other phases of gas, which yields an effective $\alpha$ less than 2.

Assuming the gravitational potential is an isothermal sphere, we have
$M(r)=2\sigma_{\rm cl}^{2}\,r/G$, where $\sigma_{\rm cl}$ is the velocity
dispersion of the clumps. Equation 2 can then be analytically solved as:

$\displaystyle v(r)=\sqrt{4\,\sigma_{\rm cl}^{2}\,{\rm ln}\Big{(}\frac{r_{\rm
min}}{r}\Big{)}+v_{\rm cl,\,\infty}^{2}\Big{(}1-\Big{(}\frac{r}{r_{\rm
min}}\Big{)}^{1-\alpha}\Big{)}}$ (3)

where $r_{\rm min}$ is the inner cutoff (or “launching”) radius that satisfies
$v(r_{\rm min})=0$, and $v_{\rm cl,\,\infty}$ = $\sqrt{2Ar_{\rm
min}^{1-\alpha}/(\alpha-1)}$ is the asymptotic maximum outflow velocity if
there were no gravitational deceleration. Note that in general the actual
$v(r)$ does not reach $v_{\rm cl,\,\infty}$ due to the gravitational
deceleration term; even the maximum radial $v(r)$ is usually several hundred
km s-1 smaller than $v_{\rm cl,\,\infty}$. Following Dijkstra & Kramer (2012),
we have fixed $\alpha$ = 1.4 and left $\sigma_{\rm cl}$ and $v_{\rm
cl,\,\infty}$ as the free parameters in this model. We set $r_{\rm min}$ , so
that $\frac{r}{r_{\rm min}}\in[1,100]$. The model is intrinsically rescalable
(i.e. increasing the size of every component in the model by any factor with
all column densities unchanged would yield an identical model) and constrains
only the ratio $\frac{r}{r_{\rm min}}$, so the following analysis applies to
$\rm Ly\alpha$ halos of varying physical sizes.

For each multiphase, clumpy model on the grid, MCRT has been performed on 104
$\rm Ly\alpha$ photon packages emitted at the center of the simulation sphere
in the form of a normalized Gaussian intrinsic spectrum
$\mathcal{N}(v,\mu=0,\sigma=\sigma_{\rm i,\,cl})$, where $\sigma_{\rm i,\,cl}$
= 12.85 km s-1 is the canonical thermal velocity dispersion of the H I gas in
the clumps at $T$ = 104 K.101010In the multiphase, clumpy model, the width of
the intrinsic spectrum is always assumed to be small and the clump velocity
dispersion is responsible for broadening the spectrum. Such a choice has the
advantage of avoiding obtaining unphysically large intrinsic line widths from
fitting the spectrum (e.g. using the shell models, see Li & Gronke 2022). The
H I clumps with a constant column density $N_{\rm HI,\,cl}$ are placed
uniformly radially, so that their number density $n_{\rm cl}\propto r^{-2}$
(i.e. mass conserving if the radial outflow velocity is constant).

Each model on the grid is further used to generate three spatially binned $\rm
Ly\alpha$ profiles by separating all the photons into three spatial bins
according to their last-scattering impact parameters: $b/b_{\rm max}\in$ (0,
0.25], (0.25, 0.50] and (0.50, 0.75], where $b_{\rm max}$ is the largest
impact parameter of the scattered $\rm Ly\alpha$ photons (see Figure 5 of Li
et al. (2022) for an illustrative schematic), and the impact parameter $b$ is
measured orthogonal to the direction of the photon’s escape trajectory. The
difference between $b_{\rm max}$ and the halo radius $r_{\rm h}$ is
negligible, and we simply

We only include the photons within 75% of $b_{\rm max}$ (or equivalently,
within the inner $\sim$ 56% of the total area) in our fitting in order to
ensure a direct comparison between the model and the data, because the S/N
$>2$ regions of the halos used for the spectra (see Section 5.1) contain on
average 58% of the total halo area. The spectra to be modeled are constructed
in the same way as the annular spectra described in Section 5.1, except that
the spaxels are divided into three radial bins with $0<r\leq 0.33r_{\rm max}$,
$0.33r_{\rm max}<r\leq 0.67r_{\rm max}$, and $0.67r_{\rm max}<r\leq r_{\rm
max}$, where $r_{\rm max}$ is the radius of the most distant spaxel in the
modeled area. When we present our modeling results later in §6.3, we consider
only the photons included in the modeling and renormalize the halo to 0.75
$b_{\rm max}$, so that $b/b_{\rm max}\in(0,\frac{1}{3}]$,
$(\frac{1}{3},\frac{2}{3}]$ and $(\frac{2}{3},1]$.

Our fitting pipeline employs the python nested sampling package dynesty
(Skilling 2004, 2006; Speagle 2020). At each visited point of the parameter
space, the pipeline executes the following three steps:

(1) calculate three binned $\rm Ly\alpha$ model spectra via linear flux
interpolation on the model grid (to circumvent doing computationally expensive
RT “on the fly”), where the flux density of the model spectrum at each
wavelength is calculated by a parameter-weighted multidimensional linear
interpolation111111Such an interpolation is carried out based on the distance
between the visited point in the parameter space and its adjacent points on
the grid (realized by the PYTHON function scipy.interpolate.interpn). of the
flux densities of the adjacent grid model spectra at the corresponding
wavelength. The three binned $\rm Ly\alpha$ model spectra are then convolved
with a Gaussian function with $\sigma=65\,\rm km\,s^{-1}$ (the KCWI line
spread function [LSF]) to mimic the finite instrumental resolution;

(2) compare each binned model spectrum to an observed $\rm Ly\alpha$ spectrum
at the corresponding impact parameter range and calculate the likelihood;

(3) sum the likelihoods of these three binned models as the likelihood of the
current set of parameters.

Each fitting run yields a posterior probability distribution (PDF) of the
model parameters. The parameter uncertainties can be further determined as
certain quantiles (e.g. 16%–84%, or 1$\sigma$ confidence intervals) of the
samples in the marginalized PDF.

###  6.2 Metal absorption line modeling

In addition to the $\rm Ly\alpha$ profiles observed at both $b=0$ and $b>0$,
the rest-UV, low-ionization metal absorption lines observed “down-the-barrel”
(i.e., at $b=0$) also encode rich information on the physical properties of
the cool gas. These metal absorption line profiles are typically “sawtooth”
shaped (e.g. Weiner et al. 2009), where the part blueward of the absorption
trough (the location of the minimum flux density) gradually decreases with
velocity while the part redward of the absorption trough increases with
velocity relatively rapidly. The blueshifted absorption at negative velocities
is produced by gas clumps with radially varying outflow velocities along the
line-of-sight, whereas the red part is mainly produced by a group of non-
outflowing, randomly moving clumps. In this work, we focus on modeling the
portion blueward of the absorption trough for the average line profile of the
Si II $\lambda$1260 and C II $\lambda$1334 transitions,121212We did not fit Si
II $\lambda$1260 and C II $\lambda$1334 separately as many of the individual
lines have fairly low S/N ratios. as we are most interested in constraining
the clump outflow kinematics. Our model is similar to the kinematic model used
by Steidel et al. (2010), but with a different clump radial velocity profile.

In our model, we first assume that the clump radial outflow velocity is
described by the same model we use for the $\rm Ly\alpha$ emission, i.e.
Equations 2 and 3 with two free parameters: the clump velocity dispersion
$\sigma_{\rm cl}$ and the asymptotic maximum clump outflow velocity $v_{\rm
cl,\,\infty}$. Assuming that the absorption lines are saturated131313The
assumption of saturation comes from the fact that in our sample, Si II
$\lambda$1260 and Si II $\lambda$1526 have similar equivalent widths (see e.g.
footnote 27 of Steidel et al. 2018). (i.e. the column densities of the
absorbing gas are so high that the depth of absorption simply reflects the gas
covering fraction), the down-the-barrel absorption line profile $I(v)$ (the
normalized, residual flux density as a function of velocity) is simply given
by

$I(v)=1-f_{c}(v)$ (4)

where $f_{c}(v)$ is the (clumpy) gas geometric covering fraction as a function
of velocity, which is the fraction of the total lines of sight of the rest-UV
emission that are intercepted by the absorbing gas. We further assume that the
gas covering fraction decreases as a function of radius, in the form of a
power law:

$f_{c}(r)=f_{c,\rm max}\left(\frac{r}{r_{\rm min}}\right)^{-\gamma}$ (5)

where $r_{\rm min}$ is the launching radius and $f_{c,\rm max}$ is the maximum
gas covering fraction that corresponds to the deepest part of the absorption
trough. $f_{c}(r)$ can then be translated into $f_{c}(v)$ using the $v(r)$
dictated by Equation 3. Note that the gas geometric covering fraction in the
$\rm Ly\alpha$ RT models, which is a function of the number density and the
physical size of the clumps (both of which may vary as a function of velocity
or radius; see Equation 2 in Dijkstra & Kramer 2012), may not be fully
consistent with the power law $f_{c}(r)$ assumed here. One may match them by
using clumps with radially-varying sizes in the RT model; we plan to explore
this option in future work.

To be consistent with the $\rm Ly\alpha$ modeling in §6.1, we fix $\alpha$ =
1.4 in the clump radial velocity profile and set $r_{\rm min}$ = 0.1 kpc with
$\frac{r}{r_{\rm min}}\in[1,100]$. Note again that only the ratio
$\frac{r}{r_{\rm min}}$ (rather than $r$ or $r_{\rm min}$ individually) is
constrained by the absorption line modeling. We then fit the observed
absorption line profiles with dynesty to determine the PDF of the four
parameters in this model: $\sigma_{\rm cl}$, $v_{\rm cl,\,\infty}$,
$f_{c,\,\rm max}$ and $\gamma$. We use flat priors for the fitted parameters
(which can vary continuously): $\sigma_{\rm cl}\in[0,120]$ km s-1, $v_{\rm
cl,\,\infty}\in[100,1500]$ km s-1, $f_{c,\,\rm max}\in[0,1]$, and
$\gamma\in[0.1,2.0]$. We restrict $\gamma$ to be no larger than 2, as
otherwise it suggests that the clumps are destroyed rapidly as they move
outwards, contradictory to the observation of metal absorption at large impact
parameters ($b\sim 100$ kpc, see Figure 21 of Steidel et al. 2010 and Rudie et
al. 2019).

Figure 12: Modeling results of the annular-averaged, spatially resolved $\rm
Ly\alpha$ spectra and of the average line profile of Si II $\lambda$1260 and C
II $\lambda$1334 observed down the barrel for Q0207-BX144 (see Appendix A for
the rest of the sample). The top row shows the best-fit models (red) to the
spatially resolved $\rm Ly\alpha$ spectra (black, with 1-$\sigma$
uncertainties shown in grey) from the inner to the outer halo. In each
subpanel of the top row, the vertical and horizontal black dashed lines
indicate the systemic redshift (determined from nebular emission lines) and
zero flux density, respectively. The middle row and the first panel of the
bottom row show a comparison between the radial trends of peak separation,
blue-to-red flux ratio, trough flux fraction, and normalized SB versus the
normalized impact parameter predicted by the best-fit models (red squares) and
measured from observation (black points, with 1-$\sigma$ uncertainties). Note
that the impact parameters may be slightly different for the model and the
data: the models are binned consistently as $b/b_{\rm max}\in(0,\frac{1}{3}]$,
$(\frac{1}{3},\frac{2}{3}]$ and $(\frac{2}{3},1]$, and while the data are
binned in the same way, the halos are asymmetric with the result that the
median distance to the spaxels included in each bin varies from object to
object. The rest of the bottom row shows the best-fit models (red) to the
average line profile (black, with 1-$\sigma$ uncertainties shown in grey) of
Si II $\lambda$1260 (blue) and C II $\lambda$1334 (orange) profiles, as well
as a comparison of clump radial outflow velocity profiles inferred from $\rm
Ly\alpha$ RT modeling (red) and metal absorption line fitting (blue hatched
patch). The shaded regions represent the velocity ranges spanned by 50 points
in the parameter space after convergence has been achieved for the fitting.
Table 5: Best-fit parameters from modeling $\rm Ly\alpha$ emission and the
rest-UV low ionization metal absorption lines

|  Best-fit Parameters ($\rm Ly\alpha$) | Best-fit Parameters (Absorption)
---|---|---
ID | $F_{\rm V}$ | ${\rm log}\,N_{\rm HI,\,{\rm cl}}$ | $\sigma_{\rm cl}$ | $v_{\rm cl,\,\infty}$ | $v_{\rm cl,\,max}$ | ${\rm log}\,n_{\rm HI,\,{\rm ICM}}$ | $\Delta v$ | $\sigma_{\rm cl,\,abs}$ | $v_{\rm cl,\,\infty,\,abs}$ | $v_{\rm cl,\,max,\,abs}$ | $f_{c,\,\rm max}$ | $\gamma$
|  |  (cm${}^{-2})$ |  (km s${}^{-1})$ |  (km s${}^{-1})$ |  (km s${}^{-1})$ |  (cm${}^{-3})$ |  (km s${}^{-1})$ |  (km s${}^{-1})$ |  (km s${}^{-1})$ |  (km s${}^{-1})$ |  |
(1) | (2) | (3) | (4) | (5) | (6) | (7) | (8) | (9) | (10) | (11) | (12) | (13)
Q0142-BX165 |  0.06${}^{+0.01}_{-0.01}$ |  18.3${}^{+0.1}_{-0.1}$ |  142${}^{+4}_{-6}$ |  748${}^{+24}_{-25}$ |  392${}^{+26}_{-27}$ |  -7.13${}^{+0.05}_{-0.07}$ |  9${}^{+4}_{-4}$ |  58${}^{+36}_{-33}$ |  835${}^{+69}_{-69}$ |  707${}^{+63}_{-46}$ |  0.8${}^{+0.1}_{-0.1}$ |  1.7${}^{+0.2}_{-0.3}$
Q0142-BX186 |  0.12${}^{+0.02}_{-0.02}$ |  18.8${}^{+0.1}_{-0.1}$ |  117${}^{+14}_{-4}$ |  533${}^{+48}_{-23}$ |  235${}^{+17}_{-4}$ |  -7.28${}^{+0.16}_{-0.13}$ |  55${}^{+10}_{-10}$ |  72${}^{+32}_{-41}$ |  429${}^{+400}_{-205}$ |  423${}^{+432}_{-157}$ |  0.6${}^{+0.3}_{-0.3}$ |  1.2${}^{+0.6}_{-0.7}$
Q0207-BX87 |  0.08${}^{+0.01}_{-0.00}$ |  18.0${}^{+0.0}_{-0.0}$ |  98${}^{+0}_{-1}$ |  617${}^{+11}_{-10}$ |  389${}^{+13}_{-12}$ |  -7.14${}^{+0.02}_{-0.02}$ |  -4${}^{+5}_{-4}$ |  61${}^{+40}_{-34}$ |  1029${}^{+311}_{-264}$ |  891${}^{+292}_{-262}$ |  0.2${}^{+0.1}_{-0.1}$ |  1.1${}^{+0.6}_{-0.6}$
Q0207-BX144 |  0.08${}^{+0.01}_{-0.01}$ |  18.4${}^{+0.1}_{-0.1}$ |  106${}^{+17}_{-21}$ |  842${}^{+40}_{-47}$ |  622${}^{+17}_{-17}$ |  -6.69${}^{+0.04}_{-0.04}$ |  -47${}^{+3}_{-3}$ |  53${}^{+37}_{-29}$ |  807${}^{+79}_{-76}$ |  684${}^{+72}_{-50}$ |  0.8${}^{+0.1}_{-0.1}$ |  1.7${}^{+0.2}_{-0.3}$
Q0449-BX110 |  0.09${}^{+0.02}_{-0.01}$ |  18.3${}^{+0.1}_{-0.1}$ |  9${}^{+11}_{-6}$ |  892${}^{+5}_{-9}$ |  815${}^{+5}_{-9}$ |  -6.89${}^{+0.04}_{-0.04}$ |  33${}^{+5}_{-5}$ |  62${}^{+36}_{-37}$ |  1051${}^{+130}_{-102}$ |  913${}^{+121}_{-84}$ |  0.8${}^{+0.1}_{-0.1}$ |  1.6${}^{+0.3}_{-0.4}$
Q0449-BX115 |  0.13${}^{+0.02}_{-0.02}$ |  18.1${}^{+0.1}_{-0.1}$ |  18${}^{+20}_{-13}$ |  550${}^{+19}_{-14}$ |  495${}^{+9}_{-9}$ |  -6.93${}^{+0.04}_{-0.04}$ |  1${}^{+3}_{-3}$ |  58${}^{+39}_{-35}$ |  636${}^{+307}_{-155}$ |  502${}^{+319}_{-131}$ |  0.7${}^{+0.2}_{-0.2}$ |  1.3${}^{+0.5}_{-0.6}$
Q0821-MD36 |  0.15${}^{+0.00}_{-0.01}$ |  18.1${}^{+0.0}_{-0.0}$ |  60${}^{+4}_{-3}$ |  506${}^{+10}_{-4}$ |  386${}^{+9}_{-10}$ |  -7.47${}^{+0.03}_{-0.02}$ |  6${}^{+4}_{-7}$ |  61${}^{+36}_{-37}$ |  870${}^{+383}_{-269}$ |  735${}^{+373}_{-249}$ |  0.3${}^{+0.3}_{-0.2}$ |  1.2${}^{+0.5}_{-0.7}$
Q1549-BX102 |  0.07${}^{+0.01}_{-0.01}$ |  17.7${}^{+0.1}_{-0.1}$ |  134${}^{+6}_{-11}$ |  882${}^{+11}_{-16}$ |  578${}^{+24}_{-19}$ |  -6.67${}^{+0.03}_{-0.04}$ |  -15${}^{+4}_{-4}$ |  58${}^{+40}_{-33}$ |  633${}^{+170}_{-125}$ |  497${}^{+166}_{-100}$ |  0.7${}^{+0.1}_{-0.2}$ |  1.4${}^{+0.4}_{-0.6}$
Q1700-BX729 |  0.16${}^{+0.00}_{-0.00}$ |  18.5${}^{+0.0}_{-0.0}$ |  123${}^{+1}_{-3}$ |  606${}^{+13}_{-6}$ |  287${}^{+16}_{-11}$ |  -6.60${}^{+0.07}_{-0.09}$ |  67${}^{+7}_{-5}$ |  44${}^{+35}_{-24}$ |  563${}^{+74}_{-65}$ |  468${}^{+44}_{-60}$ |  0.9${}^{+0.1}_{-0.1}$ |  1.7${}^{+0.2}_{-0.3}$
Q2206-BX151 |  0.10${}^{+0.02}_{-0.02}$ |  17.6${}^{+0.1}_{-0.2}$ |  142${}^{+5}_{-8}$ |  864${}^{+17}_{-23}$ |  530${}^{+14}_{-15}$ |  -6.73${}^{+0.03}_{-0.03}$ |  -73${}^{+3}_{-3}$ |  57${}^{+39}_{-33}$ |  743${}^{+258}_{-137}$ |  612${}^{+251}_{-127}$ |  0.5${}^{+0.1}_{-0.1}$ |  1.4${}^{+0.4}_{-0.6}$
Q2343-BX418 |  0.07${}^{+0.01}_{-0.01}$ |  18.0${}^{+0.0}_{-0.1}$ |  6${}^{+14}_{-4}$ |  606${}^{+19}_{-14}$ |  553${}^{+11}_{-12}$ |  -6.99${}^{+0.03}_{-0.02}$ |  34${}^{+2}_{-2}$ |  59${}^{+39}_{-35}$ |  857${}^{+165}_{-107}$ |  725${}^{+154}_{-96}$ |  0.6${}^{+0.2}_{-0.2}$ |  1.4${}^{+0.4}_{-0.5}$
Q2343-BX660 |  0.07${}^{+0.01}_{-0.01}$ |  18.4${}^{+0.1}_{-0.1}$ |  134${}^{+8}_{-19}$ |  876${}^{+17}_{-47}$ |  571${}^{+21}_{-20}$ |  -6.62${}^{+0.04}_{-0.05}$ |  -112${}^{+5}_{-5}$ |  55${}^{+37}_{-32}$ |  697${}^{+78}_{-78}$ |  569${}^{+69}_{-44}$ |  0.9${}^{+0.1}_{-0.1}$ |  1.5${}^{+0.3}_{-0.3}$

††footnotetext: Notes. Best-fit parameters (averages and 16% – 84% quantiles,
i.e., 1$\sigma$ confidence intervals) from the $\rm Ly\alpha$ and low
ionization metal absorption line (the average of Si II $\lambda$1260 and C II
$\lambda$1334) modeling. The columns are: (1) the object ID; (2) the clump
volume filling factor; (3) the clump H I column density; (4) the clump
velocity dispersion; (5) the clump asymptotic outflow velocity; (6) the actual
maximum clump radial outflow velocity; (7) the residual H I number density of
the ICM; (8) the velocity shift relative to the systemic redshift of the
source. (9) - (13) are determined from the average metal absorption line
profile. (9) the clump velocity dispersion; (10) the clump asymptotic outflow
velocity; (11) the actual maximum clump radial outflow velocity; (12) the
maximum clump covering fraction; (13) the power-law index of the clump
covering fraction function.

###  6.3 Modeling results & interpretation

Our modeling of both the spatially resolved $\rm Ly\alpha$ emission and the UV
absorption lines has achieved the following principal results:

(1) reproducing the radially varying, spatially resolved $\rm Ly\alpha$
profiles;

(2) reproducing the radial trends of several important physical quantities of
the $\rm Ly\alpha$ profiles, including the peak separation, peak flux ratio,
trough flux fraction, and SB versus the impact parameter;

(3) reconciling the clump outflow velocities inferred from $\rm Ly\alpha$
emission and metal absorption lines.

We present the modeling results for the spatially resolved $\rm Ly\alpha$
spectra and the average line profile of Si II $\lambda$1260 and C II
$\lambda$1334 for our sample in Figure 12 (using Q0207-BX144 as an example)
and Appendix A. In each panel, the top row shows the best-fit RT models (red)
to the spatially resolved $\rm Ly\alpha$ spectra (black); the middle row and
the first panel of the bottom row show a comparison between the radial trends
of peak separation, peak flux ratio, trough flux fraction, and SB predicted by
the best-fit models and measured from observations; and the rest of the bottom
row shows the best-fit models (red) to the average metal absorption line
profile (black), as well as a comparison of clump radial outflow velocity
profiles inferred from $\rm Ly\alpha$ emission and the average metal
absorption line. The best-fit parameters are summarized in Table 5, . In
Section 6.3.1 below we describe the relationships between impact parameter,
the properties of the model $\rm Ly\alpha$ profiles, and the parameters of the
model, and in Section 6.3.2 we further discuss the best-fit parameters and
relationships between them. Section 6.3.3 provides a comparison of spatially
integrated vs. spatially resolved $\rm Ly\alpha$ modeling, and we discuss
caveats to the models in Section 6.3.4.

####  6.3.1 Radial trends

Figure 13: Schematic of the escape of two $\rm Ly\alpha$ photons at low and
high impact parameters in the multiphase, clumpy model. The large circle
represents the boundary of the simulated spherical region, divided into shaded
green, red hatched, and blue regions indicating the three ranges of impact
parameters modeled. The location of the observer is indicated by the telescope
dome at the bottom, and the dotted horizontal line indicates that photons in
the blue peak arise from the near side of the halo while those in the red peak
predominantly come from the far side. The gold sun symbol represents the $\rm
Ly\alpha$ emitting source at the center, the grey clouds represent H I clumps
with random motions and radial outflows, and the small red circles represent
the diffuse, hot ICM. The impact parameters $b$ and $b^{\prime}$ are defined
as the orthogonal distance from the center to the direction of the photon
escape trajectories shown by the black solid and dashed lines. The photon that
escapes at a higher $b>b^{\prime}$ will experience several differences before
it escapes: (1) it will scatter with lower H I column densities from the
clumps, due to the decrease in the clump covering fraction at large radii; (2)
it will experience (on average) a lower projected component of the clump
outflow velocity along its traveling direction ($v_{\rm
cl,\,\parallel}<v^{\prime}_{\rm cl,\,\parallel}$, as indicated by the black
arrows near the last clump that scatters each photon); (3) it will suffer from
less absorption at line center from the ICM, due to its lower traveling
distance at the outskirts of the halo. Also note that the photon escaping at
$b^{\prime}$ passes through a clump on the near side of the halo unimpeded,
because it is out of resonance with the clump due to its previous scattering.

The modeling results show that our multiphase, clumpy model is able to
reproduce the spatially resolved $\rm Ly\alpha$ spectra fairly well,
especially for the innermost two spatial bins. In a number of cases (e.g.
Q0142-BX165, Q0207-BX87, Q0207-BX144, Q1549-BX102 and Q2343-BX660) there is a
noticeable mismatch between the model and data in the outermost bin, which may
be because the gas in the outer halo does not fully follow the outflowing
kinematics of the gas in the inner halo (e.g. due to external forces). In
general, as the impact parameter increases, the best-fit $\rm Ly\alpha$ RT
model predicts a decrease in the peak separation, an increase in the blue-to-
red peak flux ratio, and an increase in the trough flux fraction. These
predicted radial trends of peak separation and peak flux ratio are broadly
consistent with the observational data, although the exact values differ in
some cases. The increase in the trough flux fraction is also evident in almost
all objects, especially from a comparison between the innermost two spatial
bins.

From a $\rm Ly\alpha$ RT perspective, the peak separation, which reflects the
most likely frequencies at which the $\rm Ly\alpha$ photons escape, is
directly related to the $\rm Ly\alpha$ optical depth of the system. The
optical depth, which is the product of the $\rm Ly\alpha$ cross-
section141414Strictly speaking, the peak separation is also related to the gas
outflow velocity, since the $\rm Ly\alpha$ cross section depends on the
photons’ apparent frequencies in the gas frame. However, our tests have shown
that such an effect is minor compared to the one that the H I column density
has on peak separation. and the H I column density of the absorber, can
therefore be expressed as a function of the temperature and column density of
the absorber.151515For example, the peak separation of $\rm Ly\alpha$ photons
that escape from an opaque, static H I sphere is $\Delta v_{\rm peak}\simeq
320\Big{(}\frac{N_{\rm HI}}{10^{20}\,{\rm
cm^{-2}}}\Big{)}^{1/3}\Big{(}\frac{T_{\rm HI}}{10^{4}\,{\rm
K}}\Big{)}^{1/6}\,\rm km\,s^{-1}$ (Dijkstra 2014b). The blue-to-red peak flux
ratio, however, is negatively correlated with the H I gas outflow velocity as
seen by the $\rm Ly\alpha$ photons, as the blue photons are less likely to
escape since they appear closer to resonance in the reference frame of the
outflowing gas. Finally, as the absorption at the line center is mainly
produced by the ICM, the trough flux fraction is mostly set by the ICM column
density.

One can then imagine that the $\rm Ly\alpha$ photons that escape at large
impact parameters (i.e. the directions of their escape trajectories are almost
orthogonal to the radial direction) will experience the following differences
relative to photons from smaller impact parameters before they escape: (1)
experience lower H I column densities from the clumps, as the area covering
fraction of the clumps decreases at large radii due to the increase of the
physical volume of the halo; (2) encounter (on average) a lower projected
component of the clump outflow velocity along their traveling directions in
the portion of the outer halo that they pass through before they
escape;161616This is a purely geometrical effect; assuming the clump outflow
is nearly isotropic, at high impact parameters ($b\simeq b_{\rm max}$) the
maximum projected component of the clump outflow velocity along the traveling
direction of a photon goes as $v_{\rm
cl,\,\parallel,\,max}(r)=\sqrt{1-(b/b_{\rm max})^{2}}\,v_{\rm cl}(r)\simeq 0$
(Li et al. 2022). (3) suffer from lower absorption (or equivalently, “see” a
lower optical depth) at line center from the ICM in the outer halo, as on
average the distance a photon travels within the halo before it escapes at
large impact parameters is smaller than that at small impact
parameters.171717Considering the spherical geometry of the halo, the largest
distance that a photon can travel through without changing direction at impact
parameter $b$ is $\sim 2\sqrt{R^{2}-b^{2}}$. These three effects are
presumably responsible for the observed radial variation of the spatially
binned $\rm Ly\alpha$ profiles, and we illustrate them in Figure 13.

To test these hypotheses, we have designed several experiments and present
them in Figure 14. We first generate our fiducial model by setting ($F_{\rm
V}$, ${\rm log}\,N_{\rm HI,\,cl}$, $\sigma_{\rm cl}$, $v_{\rm cl,\,\infty}$,
${\rm log}\,n_{\rm HI,\,ICM}$, $\Delta v$) = (0.05, 18.5, 80, 500, -7.0, ).
Such a choice roughly corresponds to the median parameter values of the model
grid and proves to clearly demonstrate the radial variation of the peak
separation, peak flux ratio, and the trough flux fraction of the radially
binned $\rm Ly\alpha$ spectra. We then generate three test models for
comparison by modifying the configuration of the fiducial model in specific
ways. In Model I, we adjust the spatial distribution of the clumps: instead of
placing the clumps radially uniformly, we place more clumps at large radii so
that the number density of the clumps $n_{\rm cl}(r)\simeq$ constant. In Model
II, we change the direction of the clumps’ outflow velocity from radial to
tangential , so that the projected component of the clump outflow velocity
along the traveling direction is no longer preferentially small for photons
that escape at high impact parameters. In Model III, we increase the number
density of the ICM by a factor of 20 in the outer 60% of the halo radius in
order to offset the shorter photon traveling distance at large radii. As shown
in Figure 14, in Model I, the peak separation of the three binned $\rm
Ly\alpha$ model spectra is now roughly constant; in Model II, the significant
increase in the blue-to-red peak flux ratio is no longer present, yet a slight
decrease towards the outskirts is seen; and in Model III, the trough flux
fractions are all much closer to zero. Therefore, we conclude that these
experiments strongly support our above explanation for the radial trends of
the peak separation, peak flux ratio and trough flux fraction.

Figure 14: Experiments designed to test our hypotheses for the differences
between $\rm Ly\alpha$ photons that escape at low and high impact parameters.
In each of the four subpanels, three binned model $\rm Ly\alpha$ spectra are
shown according to their last-scattering impact parameters: $b/r_{\rm
h}\in(0,\frac{1}{3}]$ (green solid), $(\frac{1}{3},\frac{2}{3}]$ (red dash-
dotted) and $(\frac{2}{3},1]$ (blue dotted), where $r_{\rm h}$ is the radius
of the modeled halo. _Left_ : The fiducial model with ($F_{\rm V}$, ${\rm
log}\,N_{\rm HI,\,cl}$, $\sigma_{\rm cl}$, $v_{\rm cl,\,\infty}$, ${\rm
log}\,n_{\rm HI,\,ICM}$) = (0.05, 18.5, 80, 500, -7.0). _Second from left_ :
Model I, in which more clumps are placed at large radii so that the number
density of the clumps $n_{\rm cl}(r)\simeq$ constant. _Third from left_ :
Model II, in which the clump radial velocity is set to be tangential, so that
the projected component of the clump outflow velocity along the traveling
direction is no longer preferentially small for photons that escape at high
impact parameters. _Right_ : Model III, in which the number density of the ICM
is increased by a factor of 20 in the outer 60% of the halo radius. In each of
the three test models, the change in the model configuration offsets the
corresponding spatial variation of the $\rm Ly\alpha$ spectral morphology
(i.e. peak separation, peak flux ratio and trough flux fraction), hence
supporting our explanation.

Incidentally, our model has also reproduced the decreasing trend of $\rm
Ly\alpha$ SB versus impact parameter, with only a few exceptions (e.g.
Q0142-BX186 and Q1700-BX729). These two objects, which have a more gradual
decline in SB, are the faintest objects in the sample, with the smallest
fraction of the total halo area used for the spatially resolved $\rm Ly\alpha$
modeling. This overall consistency adds further credence to our multiphase,
clumpy RT model.

####  6.3.2 Best-fit parameters

One of the most interesting discoveries from our modeling is that the clump
outflow velocities inferred from $\rm Ly\alpha$ emission and the low
ionization metal absorption lines can be mutually consistent, with typical
values of $\sim 400-600$ km s-1 obtained for both . The mismatch between the
gas outflow velocities inferred from $\rm Ly\alpha$ and from metal absorption
lines has been a long-standing problem. For example, it is reported that the
$\lesssim 150$ km s-1 outflow velocities of the shell model required to match
the $\rm Ly\alpha$ profiles of local starburst and green pea galaxies are much
lower than the $\gtrsim 300$ km s-1 characteristic velocities of the metal
absorption lines (e.g. Leitherer et al. 2013; Orlitová et al. 2018). The high
outflow velocity regime of $\rm Ly\alpha$ RT models has been little explored,
possibly due to the belief that the $\rm Ly\alpha$ photons will be seen as out
of resonance by the fast moving gas and will therefore not scatter (e.g.
Verhamme et al. 2015). However, we observe an interesting pattern in our
multiphase, clumpy model: for a typical double-peaked $\rm Ly\alpha$ profile,
as the clump outflow velocity increases, the blue-to-red peak flux ratio (or
the “level of symmetry”) first decreases and then increases, until the clump
outflow velocity is so large that all the photons are completely shifted out
of resonance as seen by the gas. This pattern, as shown by an example in
Figure 15, suggests the possibility of matching the observed asymmetric $\rm
Ly\alpha$ profiles in the high outflow velocity regime ($v_{\rm
cl,\,max}\gtrsim 400-600\rm\;km\,s^{-1}$).

In our sample, consistency (accounting for uncertainties) between the clump
outflow velocities inferred from $\rm Ly\alpha$ and metal absorption lines is
achieved in 8 / 12 objects.181818We define two velocity profiles as being
consistent if they have a non-negligible overlap at $r>r_{\rm min}$. Such a
high success rate demonstrates the feasibility of matching both the observed
$\rm Ly\alpha$ and metal absorption line profiles simultaneously with one
clump radial velocity profile. Among the four inconsistent cases, two
(Q0207-BX87 and Q2343-BX418) have relatively irregular and noisy absorption
line profiles that yield a broad range of velocities, whereas in the other two
cases (Q0142-BX165 and Q1700-BX729), . We also note that exact matches between
the $\rm Ly\alpha$ and absorption-line-inferred outflow velocities are not
necessarily expected because the transitions probe somewhat different gas: the
absorption lines are purely a line-of-sight measurement that probes the gas
only on the near side of the halo, while the $\rm Ly\alpha$ results
incorporate gas on the far side and at large impact parameters that is not
seen in absorption.

---
Figure 15: Examples of $\rm Ly\alpha$ model spectra with different clump
outflow velocities showing the pattern in the change of the blue-to-red peak
flux ratio. Five models with ($F_{\rm V},{\rm log}\,N_{\rm HI,\,cl},{\rm
log}\,n_{\rm ICM},\sigma_{\rm cl})=(0.04,17.5,-6.5,125)$ and $v_{\rm
cl,\,\infty}=(500,600,700,800,900)$ are shown with different colors and
linestyles. _Upper_ : spatially integrated $\rm Ly\alpha$ model spectra with
different $v_{\rm cl,\,\infty}$ values. As $v_{\rm cl,\,\infty}$ increases,
the average clump radial outflow velocity increases, and the blue-to-red peak
flux ratio first decreases (comparing the black and green curves) and then
increases (comparing the red, blue and orange curves). _Lower_ : The
corresponding clump radial velocity profiles for different $v_{\rm
cl,\,\infty}$. Note that $v_{\rm cl}(r)$ and $v_{\rm cl,\,\infty}$ are
positively correlated, but $v_{\rm cl}(r)$ is always smaller than (typically
by several hundred km s-1) $v_{\rm cl,\,\infty}$ due to the effect of
gravitational deceleration.

The best-fit radial velocity profile of the clumps in the multiphase, clumpy
model typically exhibits a rapid acceleration phase to $v_{\rm cl}=v_{\rm
cl,\,max}$ within $1\lesssim\frac{r}{r_{\rm min}}\lesssim 10$ followed by a
gradual deceleration191919Note, however, that such a deceleration phase is not
preferred by the absorption line modeling as it will break the one-to-one
relation between $r$ and $v_{\rm cl}$ and yield a pathological absorption line
profile $I(v)$. (or $v_{\rm cl}\simeq$ constant) phase at $\frac{r}{r_{\rm
min}}\gtrsim 10$. The decline in the outflow velocity and possible transition
to an inflow are physically expected due to the increasing importance of
gravitational deceleration at large radii, and have been explored in previous
works (e.g. Chen et al. 2020); however, the exact location of the transition
is model-dependent and may need additional observational constraints.

We also note that there is a significant velocity difference between the
outflowing cool clumps and the static hot ICM in the best-fit models, which is
at odds with the traditional “hot wind entrains (and co-outflows with) the
cold gas” paradigm (see e.g. Gronke & Oh 2018, 2020, and references therein).
It is possible, however, that the interaction between the hot phase and the
$\rm Ly\alpha$ photons is dominated by the decelerated, semi-static hot gas,
as suggested by the deep troughs at line center in the observed $\rm Ly\alpha$
profiles. A larger sample with more diverse $\rm Ly\alpha$ morphologies will
be helpful in assessing the impact of an outflowing hot gas component in the
future.

We next turn to the other best-fit parameters of the models. For the $\rm
Ly\alpha$ modeling, the best-fit clump volume filling factors ($F_{\rm V}$)
range from 0.06 to 0.16 (corresponding to $\sim 5-10$ clumps on average per
line-of-sight202020Note that the number of clumps per line-of-sight and the
associated gas covering fraction both decrease with $r$ due to the increase of
the halo volume at large $r$ (cf. Figures 15 and 16 of Rudie et al. 2012).),
and the best-fit clump column densities ($N_{\rm HI,\,cl}$) range from $\sim
10^{17.6}$ to $10^{18.8}\,\rm cm^{-2}$. The total H I column densities
($N_{\rm HI,\,total}\simeq\frac{4}{3}f_{\rm cl}N_{\rm HI,\,cl}=(r_{\rm
h}/r_{\rm cl})\,F_{\rm V}N_{\rm HI,\,cl}$, Gronke et al. 2016) of the best-fit
models range from $\sim 10^{18.5}$ to $10^{19.9}\,\rm cm^{-2}$. Here $N_{\rm
HI,\,total}$ represents the inferred total H I column density of the modeled
halo that a $\rm Ly\alpha$ photon typically interacts with, either via
scattering or free-streaming; the scattered, out-of-resonance $\rm Ly\alpha$
photons may stream through the high-velocity, outflowing clumps without
scattering (Gronke et al. 2017).

The residual H I column densities of the hot, diffuse ICM ($N_{\rm
HI,\,ICM}\simeq n_{\rm HI,\,ICM}\,r_{\rm h}$) range from $\sim 10^{15}$ to
$10^{16}\,\rm cm^{-2}$. Such column densities are much smaller than those
within the clumps, but are necessary to produce the absorption trough at line
center, and may serve as optically-thin channels for LyC escape along lines of
sight that have relatively few H I clumps. The best-fit systemic redshifts of
the $\rm Ly\alpha$ sources are mostly consistent with the systemic redshifts
determined from nebular emission lines ($|\Delta v|<50\;\rm km\,s^{-1}$ for 8
/ 12 objects). The best-fit clump velocity dispersions ($\sigma_{\rm cl}$) are
all smaller than $150\;\rm km\,s^{-1}$ and span a similar range to the
observed nebular emission line widths ($\sim 50-120\;\rm km\,s^{-1}$). We
compared the best-fit $\sigma_{\rm cl}$ values with the MOSFIRE $H$-band ([O
III] and $\rm H\beta$) and $K$-band ($\rm H\alpha$) nebular emission line
widths (corrected for instrumental LSF), but did not find any significant
correlation.

For the metal absorption line modeling, the best-fit clump velocity
dispersions212121Note that here the velocity dispersions are determined
independently from the $\rm Ly\alpha$ modeling. In fact, they are not very
well constrained (i.e. flat posterior) by the absorption line data, as the
acceleration term is preferred to be dominant. are all smaller than $75\;\rm
km\,s^{-1}$, suggesting that the gravitational deceleration only plays a minor
role compared to the acceleration forces. The clump outflow velocities are
high, mostly $\gtrsim 500\rm\;km\,s^{-1}$, and generally correspond to the
velocity where the blue side of the absorption line profile meets the
continuum. The maximum clump covering fractions $f_{\rm c,max}$ range from
$\sim$ 0.2 to 0.9, depending on the minimum flux density of the absorption
line profile. The power-law indices of the clump covering fraction function
($\gamma$) range from $\sim$ 1.1 to 1.7, corresponding to a mass-conserving-
like (or more gradual) decrease in the number density of the clumps. A
$\gamma$ smaller than 2 may suggest that the clumps expand as they move
outwards (e.g. due to the decrease of thermal/radiation pressure at large
radii), because if the clumps are uniformly distributed radially and their
sizes remain constant at different radii, $\gamma$ will be exactly 2 due to
the geometric volume increase at large radii ($dV\propto 4\pi r^{2}dr$).

We have also checked if any correlations exist between the best-fit clump
outflow velocities and the host galaxy properties such as stellar mass and
SFR, as these are expected to be correlated due to the causal relation between
stellar feedback and galactic outflows (e.g. Martin 2005; Rupke et al. 2005;
Weiner et al. 2009; Chen et al. 2010; Martin et al. 2012; Rubin et al. 2014;
Chisholm et al. 2015; Heckman et al. 2015; Trainor et al. 2015). Specifically,
we tested for correlations between the _actual_ maximum clump radial
velocities $v_{\rm cl,\,max}$ inferred from the $\rm Ly\alpha$ and absorption
line modeling versus the stellar masses, SFRs and sSFRs of the host galaxies.
We find that all three correlations are insignificant and have considerable
scatter. Such a null result is unsurprising, however, as our sample is
intentionally restricted to low-mass galaxies with high SFR and sSFR values
and therefore has a limited dynamic range by design. We will revisit these
correlations with larger and more well-rounded samples in future work.

####  6.3.3 Advantages of spatially resolved Ly$\alpha$ modeling

In this section, we demonstrate the advantages of spatially resolved $\rm
Ly\alpha$ modeling by comparing it to the spatially integrated $\rm Ly\alpha$
modeling that has typically been carried out in previous works. Assuming that
for a $\rm Ly\alpha$-emitting source of interest, only a spatially integrated
$\rm Ly\alpha$ spectrum within a certain aperture can be obtained (e.g. due to
the unavailability of IFU observations), we consider the following two
scenarios: (1) the spatially integrated spectrum corresponds to the $\rm
Ly\alpha$ emission from only the central region, typical of observations using
a slit or other small aperture (the spectra we extracted in Section 3 and
showed in Figure 1 belong to this category); (2) the spatially integrated
spectrum is extracted from a larger aperture that also includes the $\rm
Ly\alpha$ emission from a significant portion of the extended halo. For
exploratory purposes, we model scenario (1) with spatially integrated
multiphase, clumpy models in which all the emitted photons are included in the
emergent spectra, assuming that we are completely unaware of any spatial
variation of the $\rm Ly\alpha$ emission. We model scenario (2) with spatially
integrated models that include the photons with $b/r_{\rm h}\lesssim 75\%$,
assuming that we are aware that the data only represent part of the extended
halo and should be compared to a corresponding fraction of the modeled halo.
This is equivalent to merging the 3-bin spectra for both the data and the
models in the spatially resolved modeling routine that we described in Section
6.1.

For scenario (1), we find that the best-fit clump outflow kinematics (namely
the $\sigma_{\rm cl}$ and $v_{\rm cl,\,\infty}$ values) are similar in both
the spatially integrated and resolved modeling, but the required clump volume
filling factors (and hence the covering factor) and ICM column densities are
higher, on average, in the spatially integrated $\rm Ly\alpha$ modeling. This
is mainly because in the observed spatially integrated spectra, the trough
depth at line center is similar to that of the innermost binned spectra used
in the spatially resolved modeling, as they correspond to similar regions of
the halos. In contrast, the trough depth at the line center of a spatially
integrated $\rm Ly\alpha$ model spectrum lies between that of its
corresponding innermost and outermost binned model spectra due to the radial
variation of the profile (see Section 6.3.1). Therefore, larger clump volume
filling factors (which contribute to the total H I column densities) and ICM
column densities are required to reproduce the deep troughs in the spatially
integrated $\rm Ly\alpha$ profiles.

A quantitative comparison of the best-fit total H I column densities from the
clumps and the H I densities in the ICM for the spatially resolved and
scenario (1) models is shown in Figure 16, with the darker and fainter points
indicating the resolved and spatially integrated models from scenario (1)
respectively. We plot the total $N_{\rm HI}$ and $n_{\rm HI,\,ICM}$ versus
properties measured from the integrated spectra, and discuss the comparison
further in Section 7 below. We find that values of total $N_{\rm HI}$ from
spatially integrated modeling of the central region are larger on average by a
factor of 1.5, while $n_{\rm HI,\,ICM}$ is larger by at least a factor of 1.9,
and likely significantly more because more than half of the sample requires
values of $n_{\rm HI,\,ICM}$ higher than the maximum value allowed by the
model grid. The overestimation of $n_{\rm HI,\,ICM}$ in the spatially
integrated models manifests as an overestimation of the depth of the trough
between the peaks, which is due to the omission of spatial information on the
outer halo.

For scenario (2), we find that the best-fit parameters of the spatially
resolved and integrated modeling are fully consistent with each other. This
result is probably unsurprising, as a reasonable match between all three bins
of model and data should still hold if the bins are merged for both the models
and the data. However, we stress that this result does not indicate that the
spatially resolved modeling is no longer necessary, as we would not have found
that the radial trends of peak separation, peak ratio, trough flux ratio, and
SB can all be reasonably well-matched by the same best-fit model if we had not
separated the photons into different spatial bins and modeled the $\rm
Ly\alpha$ profiles in a spatially resolved manner.

In short, our experiments in this section suggest that although spatially
integrated modeling may be used to crudely extract certain global properties
of the CGM, it tends to either lose information about the outer regions of the
halos and overestimate the neutral hydrogen content encountered by $\rm
Ly\alpha$ photons, or fail to account for the radial variation of the $\rm
Ly\alpha$ morphological properties. In comparison, spatially resolved $\rm
Ly\alpha$ modeling has the advantage of fully leveraging the spatial variation
in the $\rm Ly\alpha$ halo as observed by integral field unit spectrographs
such as KCWI and quantifying the corresponding spatial changes of the physical
parameters of the CGM. The overall good match in radial trends between the
spatially resolved data and models provides a reassuring check on the validity
of the multiphase, clumpy model.

Figure 16: Comparison of results from the radiative transfer models with
properties of the spatially integrated spectra of the central regions of the
galaxies shown in Figure 1. Top row, left to right: the $\rm Ly\alpha$ blue-
to-red flux ratio, peak separation, trough flux fraction $f_{\rm tr}$ and mean
low ionization absorption line equivalent width vs. the total H I column
density. Middle row: the same four spectral quantities vs. the residual H I
density in the ICM. Bottom row, left to right: $\rm Ly\alpha$ blue-to-red flux
ratio, peak separation, and $f_{\rm tr}$ vs. mean low ionization absorption
line equivalent width. In the top two rows the darker points show the results
of our best-fit spatially resolved modeling, while the fainter points show the
results of modeling the single, spatially integrated line profiles. The lower
corner of each panel gives the $p$-value resulting from a Spearman correlation
test using the spatially resolved models only. Values with $p<0.1$ are
highlighted in red.

####  6.3.4 Caveats

There are several important caveats to this work. First, we did not include
the effect of dust (but note that the dust extinction of our sample is
typically small), which means that all of the emitted $\rm Ly\alpha$ photons
will eventually escape from the simulation region and contribute to the
emergent model spectra. Considering that the actual $\rm Ly\alpha$ escape
fraction is always smaller than one (when it is robustly measured; see Table 3
and discussion in Section 3), we essentially assumed that the observed
frequency distribution of $\rm Ly\alpha$ photons is representative of the $\rm
Ly\alpha$ photons that escape in all directions. The validity of such an
assumption requires further scrutiny.

Second, we used spherically symmetric RT models to model the angularly
averaged $\rm Ly\alpha$ profiles of asymmetric halos, so the results should be
interpreted as average parameters within the modeled region. We have also
experimented with modeling the spatially resolved $\rm Ly\alpha$ profiles
along the directions of maximum and minimum peak ratio and peak separation
gradients (see Section 5.2), but did not find any significant dependence of
the model parameters on these higher order spatial variations. This is mainly
because the best-fit model is primarily constrained by the spectra of the two
innermost bins, which have higher S/N, whereas the spectrum of the outermost
bin may contribute strongly to the measured gradients but does not put strong
constraints on the model parameters. Development of anisotropic RT models may
shed light on this problem,

Last but not least, some of the assumptions in our models are inevitably over-
simplified. For example, we assumed a two-component model with temperatures of
104 K and 106 K, whereas in reality H I absorbers at intermediate temperatures
should exist (Rudie et al. 2019). The H I column densities and the physical
sizes of the clumps are also simplistically assumed to be constant in the
multiphase, clumpy model. Moreover, the actual motion of the clumps in the CGM
may be more complicated than the idealistic kinematic model we employed (see
e.g. Fielding & Bryan 2022). We plan to upgrade our models in future work.

##  7 Summary and discussion

We have presented KCWI integral field spectroscopy and radiative transfer
modeling of spatially extended $\rm Ly\alpha$ emission in a sample of 12
relatively low mass ($M_{\star}\sim 10^{9}$ $\rm M_{\odot}$), extreme emission
line galaxies at median redshift $z=2.3$. As described in Section 2.1, the
targets are primarily selected based on nebular emission line ratios
indicating high ionization and low metallicity, and all are previously known
$\rm Ly\alpha$-emitters. The sample galaxies have specific star formation
rates $\sim 4$ times larger than than that of their $z\sim 2$ parent sample,
and may more closely resemble galaxies at earlier epochs of cosmic history.
Our primary results are as follows:

  1. 1.

All of the galaxies show strong, double-peaked $\rm Ly\alpha$ emission (see
Section 3 and Figure 1) and spatially extended $\rm Ly\alpha$ halos, with
luminosities ranging from $3\times 10^{42}$ to $3\times 10^{43}$ erg s-1 and
radii between 16 and 30 kpc (Figure 2).

  2. 2.

We fit double asymmetric Gaussian profiles to the $\rm Ly\alpha$ emission of
individual spaxels and small Voronoi bins in each halo, as described in
Section 4 and shown in Figure 4, and measure the flux ratio of the blue and
red peaks and the peak separation velocity $\Delta v_{\rm peak}$ for each
spaxel or bin. Maps of the peak ratio and separation are presented in Figures
5 and 6. The halos show significant azimuthal variation, but the blue-to-red
flux ratio tends to increase at larger radii and regions of narrower peak
separation are usually found in the outer halo. The peak ratio and separation
are anti-correlated for half of the sample, but this is likely driven by the
underlying tendency of both to change with radius (Figures 7 and 8).

  3. 3.

We also construct spatially averaged $\rm Ly\alpha$ profiles, in order to
identify general trends and measure the profiles to larger radii. We first
construct azimuthally averaged spectra binned as a function of radius (Section
5.1 and Figure 9), and again measure the peak ratio and separation in each
annular region as well as $f_{\rm tr}$, the fraction of total flux escaping
within $\pm 100$ km s-1 of the trough between the peaks (Figure 10). The blue-
to-red flux ratio increases consistently with radius for most objects in the
sample, with a typical central value of $\sim 0.2$; all objects that can be
measured at a radius $\gtrsim 16$ kpc have peak flux ratios $>0.6$ at that
radius. $f_{\rm tr}$ also increases with radius for most of the sample. Trends
with peak separation are more complex, but the typical central peak separation
is $\sim 600$ km s-1, with a moderate decrease toward the outer halo.

  4. 4.

Because the annular binned spectra wash out the significant azimuthal
variations in the line profiles, we also construct binned spectra of 60∘
angular regions designed to maximize the gradients in peak ratio and
separation from the center to the outer halo, using seven of the brightest
galaxies in the sample (Section 5.2 and Figure 11). These spectra show that
all of the halos have sightlines for which the peak ratio increases (typically
from $\sim 0.2$ to $\sim 1$) or the peak separation decreases (typically from
$\sim 600$ – 700 to $\sim 300$ – 400 km s-1) with radius. In all cases,
however, the regions of maximum peak ratio increase and maximum peak
separation decrease do not overlap. We also construct spectra designed to
minimize the gradients in peak ratio and separation, finding that most halos
also have regions for which the changes in peak ratio and separation with
radius are relatively small.

  5. 5.

Using a new suite of $\rm Ly\alpha$ radiative transfer simulations, we model
the spatially resolved $\rm Ly\alpha$ profiles in three radial bins with
multiphase, clumpy models with radially-varying outflow velocities (Section
6.1). These models are broadly successful in reproducing the observed line
profiles, as well as the radial trends of peak flux ratio, peak separation,
and trough flux fraction (Figures 12 and 17–22). The clumps reach a typical
maximum velocity of $\sim 500$ km s-1 and have H I column densities of $\sim
10^{17.6}$ to $10^{18.8}$ cm-2, while the total $N_{\rm HI}$ of the best-fit
models ranges from $\sim 10^{18.5}$ to $10^{19.9}$ cm-2. The clumps are
embedded in a hot inter-clump medium with residual $N_{\rm HI,\,ICM}\sim
10^{15}$ – $10^{16}$ cm-2. Best-fit parameters of the models are given in
Table 5.

  6. 6.

We find that the trend in $\rm Ly\alpha$ peak separation with radius is
primarily governed by the H I column density, as photons that escape at larger
radii are able to do so with a smaller velocity shift because they experience
lower H I column densities from the clumps before they escape due to the
decrease in clump covering fraction with radius. The $\rm Ly\alpha$ peak ratio
depends on the line-of-sight velocity, with the result that the variation in
peak ratio with radius is largely a geometric effect as the projected
component of the outflow velocity along the line of sight decreases with
increasing impact parameter (Figure 13). The depth of the trough (or the
trough flux fraction, $f_{\rm tr}$) between the two peaks primarily depends on
the residual neutral H I density of the ICM. We show the results of
experiments designed to test these conclusions in Figure 14, and further
explore the relationship between outflow velocity and peak ratio in Figure 15.

  7. 7.

We self-consistently model the mean low ionization absorption line profile of
each object, employing the same radially varying velocity model used for the
$\rm Ly\alpha$ emission and a radially decreasing gas covering fraction
(Section 6.2 and Figures 12 and 17–22). Typical clump maximum outflow
velocities inferred from the absorption line profiles are $\gtrsim 500$ km
s-1, in broad agreement with the velocities inferred from $\rm Ly\alpha$;
exact matches may not be expected because the down-the-barrel UV spectra and
the radially binned $\rm Ly\alpha$ emission are not probing entirely the same
regions of the halos. This agreement alleviates a long-standing discrepancy
between outflow velocities inferred from $\rm Ly\alpha$ shell models and the
UV absorption lines.

  8. 8.

Finally, we compare the results of the spatially resolved $\rm Ly\alpha$
modeling with those obtained from applying the same model to single, spatially
integrated $\rm Ly\alpha$ profiles, using both a small aperture capturing only
the brightest region (scenario 1) and a larger aperture encompassing most of
the halo (scenario 2). We find that modeling the integrated central profile
(scenario 1) results in higher inferred values for both the total H I column
density and the neutral component of the ICM, largely because the spatially
integrated modeling does not account for the decrease in the depth of the
trough between the peaks at larger radii; this decrease in depth reflects the
lower neutral hydrogen content experienced by photons that escape from larger
radii and indicates that some photons may escape at line center in the outer
halo. The best-fit parameters obtained from modeling a larger aperture in
scenario (2) are consistent with those from the spatially resolved modeling,
but fail to capture the trends in the $\rm Ly\alpha$ profile with radius and
the physical insights these variations provide.

Our observations and modeling suggest a self-consistent physical picture of
the CGM of this sample of $z\sim 2$ star-forming galaxies: a multiphase,
clumpy medium in which cool ($\sim 10^{4}$ K), outflowing gas clumps are
embedded in a hot ($\sim 10^{6}$ K), highly ionized, diffuse medium with low-
density residual H I. The clumps typically have H I column densities of $\sim
10^{18}\,\rm cm^{-2}$ and provide a total column density of $\sim 10^{19}\,\rm
cm^{-2}$, and the $\rm Ly\alpha$ photons “solve the maze” by being resonantly
scattered by and free-streaming through the clumps until they escape. The cool
clumps also have random velocity dispersions of $\sim 100\,\rm km\,s^{-1}$,
and are accelerated to high radial outflow velocities of $\gtrsim
500\rm\,km\,s^{-1}$ at large impact parameters, which give rise to both the
asymmetric $\rm Ly\alpha$ profiles and broad low ionization metal absorption
lines. The hot ICM is nearly static and has a low total H I column density
($\sim 10^{15}$ – $10^{16}\,\rm cm^{-2}$), but is essential to shaping the
emergent double-peaked $\rm Ly\alpha$ profiles as it provides additional
scattering that produces the absorption trough at line center.

###  7.1 Central Ly$\alpha$ profiles and LyC escape

With this physical model of the CGM in mind, we revisit the spatially