The Circumgalactic Medium of Extreme Emission Line Galaxies at z ∼ 2 z ∼ 2 z ∼ 2 : Resolved Spectroscopy and Radiative Transfer Modeling of Spatially Extended Ly ααα Emission in the KBSS-KCWI Survey

The resonantly scattered Ly α 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 Ly α emission in 12 relatively low-mass ( M ⋆ ∼ 10 9 M ⊙ ) z ∼ 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 Ly α 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 Ly α 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.


INTRODUCTION
Corresponding author: Dawn K. Erb erbd@uwm.edu* Based 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.
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 Lyα transition of hydrogen, and deep observations have now revealed that the diffuse, distant universe is aglow with Lyα 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 Lyα emission surrounding typical star-forming galaxies at high redshifts came from stacked, narrowband images (Steidel et al. 2011;Momose et al. 2014Momose et al. , 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 ≲ z ≲ 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 Lyα emission.Perhaps the most straightforward of these is the resonant scattering of Lyα photons produced in galaxies by neutral hydrogen in the CGM (Zheng et al. 2011;Kusakabe et al. 2019;Byrohl et al. 2021), with the Lyα profile then reflecting the kinematics and geometry of the CGM gas.Other possible sources of the Lyα 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 Lyα 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 Lyα properties of gaseous halos are primarily influenced by galactic outflows within ∼ 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 Lyα 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 Lyα 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 Lyα 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 Lyα 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 Lyα 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 Lyα 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 Lyα 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 Lyα 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 Lyα profiles of large numbers of spatially integrated spectra, generally by modeling the outflow as a spherical, expanding shell (e.g.Verhamme et al. 2008;Verhamme et al. 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 Lyα 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, Lyα 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 Lyα profile.Such a clumpy model converges to the monolithic shell model in the limit of being "very clumpy" (i.e.having ∼ 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 KCWIobserved Lyα profiles of several regions in the z = 3.1 Lyα 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 Lyα profiles have significant residual fluxes at the line center, which correspond to relatively few clumps per lineof-sight and low residual H I density in the ICM.In addition, the very broad Lyα 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 Lyα emission.
Follow-up work by Li et al. (2022) modeled the Lyα profiles of another z = 3.1 Lyα 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 Lyα profiles.They found that the Lyα 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 lineof-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 Lyα spectra with the multiphase, clumpy model.
We analyze the spectral properties of spatially extended Lyα emission for a sample of 12 relatively low-mass, lowmetallicity galaxies at z ∼ 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 Lyα 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 Lyα emission offer the possibility of obtaining constraints on the physical conditions in the multiphase CGM gas.The double-peaked nature of Lyα emission from highly ionized sources also provides additional constraints on the models; all 12 of our targets have doublepeaked profiles, which we quantify in both individual spaxels and binned regions before modeling the results with state-ofthe-art radiative transfer codes.
We describe our sample selection, observations, and data reduction in Section 2, and measure the global properties of the Lyα emission in Section 3. In Section 4 we quantify the Lyα 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 Lyα emission and the restframe 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 , Ω m = 0.31, and Ω Λ = 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.

Sample selection and properties
Target selection for this study was motivated by the simultaneous goals of characterizing the CGM in relatively lowmass, extreme emission line galaxies and improving our ability to extract physical information from double-peaked Lyα emission.Because low-mass, low-metallicity, and highly ionized galaxies tend to exhibit strong, double-peaked Lyα 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 ∼ 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 Lyα properties of z ∼ 2 galaxies with extreme nebular emission line ratios placing them in the upper left corner of the [N II]/Hα vs. [O III]/Hβ "BPT" diagnostic diagram (Baldwin et al. 1981), with log([N II]/Hα) ≤ −1.1 and log([O III]/Hβ) ≥ 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 ∼ 2 galaxies in our sample have typical metallicities 12 + log(O/H) ≈ 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 Lyα-emitters (LAEs) in the z ∼ 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 Hα, and Q2206-BX151 has only Hα observations).The median redshift of the sample is z med = 2.32, and all of the targets fall above the canonical LAE threshold, with rest-frame equivalent width W Lyα > 20 Å measured from long-slit spectroscopy.
The sample galaxies also have very high equivalent width [O III] λ5008 emission, with W [O 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 ∼ 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 * UV = −20.70 at z ∼ 2.3 (Reddy & Steidel 2009) and median UV slope β med = −1.87.The median stellar mass of the sample is 1.5 × 10 9 M ⊙ , from modeling the spectral energy distributions with the BPASSv2.2stellar 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 M ⊙ and −1.35 between 0.1-0.5 M ⊙ .This median stellar mass implies a halo mass of ∼3 × 10 11 M ⊙ (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 ∼ 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 ) cont = 0.05.
We use the Hα/Hβ ratio and the SMC extinction law to correct the nebular emission lines for internal reddening, finding median E(B −V ) neb = 0.14.Two galaxies in the sample (Q0207-BX87 and Q2343-BX660) have Hα/Hβ 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 ) 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 ) cont = 0.0 for both Q0449-BX115 and Q2206-BX151.
Star formation rates are computed from the dust-corrected Hα luminosity using the calibration of Theios et al. (2019), who calculate the conversion between SFR and Hα luminosity for the BPASSv2.2stellar models used for the SED fitting described above.The resulting SFRs range from 2 to 25 M ⊙ yr −1 , with a median of 14 M ⊙ 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 ) cont is typically smaller than E(B−V ) neb , so for these two objects we have potentially underestimated the extinction correction and therefore also the SFR (in fact E(B − V ) 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 Hα 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 ∼ 2 slit loss correction for NIRSPEC observations of Hα emission at z ∼ 2).
From the stellar masses and SFRs we calculate the specific star formation rate, sSFR ≡ SFR/M ⋆ , 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 ∼ 2 SFR-stellar mass relation (Reddy et al. 2012;Whitaker et al. 2014, but note that sam- KBSS sample (Strom et al. 2017).High [O III]/Hβ and O32 are both associated with high equivalent width Lyα emission (e.g.Nakajima et al. 2016;Trainor et al. 2019).
We calculate the electron density from the [O II]λ3729/[O II]λ3727 ratio, which has a median value of 1.01, slightly lower than the ratios of 1.13-1.16found for composite spectra of KBSS galaxies at z ∼ 2 (Steidel et al. 2014;Strom et al. 2017).This lower ratio indicates that many of the galaxies in our sample have higher than  2. In summary, the galaxies studied here are relatively lowmass, 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 ∼ 2, but may more closely resemble galaxies observed in the reionization era (e.g.Stark et al. 2017).

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. ′′ 5 × 20.′′ 4 and spectral resolution R ≈ 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 • .

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. 2 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. ′′ 69 × 4. ′′ 6 × 100 Å.This median-filtered cube was subtracted from the original cube to remove lowfrequency 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-whitelight 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. ′′ 3 × 0. ′′ 3 × 1 Å.The resampling was conducted using the "drizzle" method in the Montage package,3 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. (3)

GLOBAL Lyα MEASUREMENTS
In this section we describe the global properties of the Lyα emission measured by KCWI.We begin with onedimensional 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. ′′ 9 to 2. ′′ 2, with all but three within 0. ′′ 3 of the sample median of 1. ′′ 6.At the median redshift of the sample, these spectra cover rest-frame wavelengths ∼1060-1660 Å, and include a number of interstellar absorption lines which we model along with the Lyα emission in Section 6.1.
We show the continuum-subtracted Lyα 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 Lyα-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 Lyα equivalent widths are given in the column labeled W 1D Lyα in Table 3.As previously known and by design, all are above the canonical Lyα-emitter threshold of W Lyα > 20 Å.
We next create pseudo-narrowband continuum-subtracted Lyα surface brightness images for each object in the sample.We first identify the spatial peak of the Lyα emission in each data cube, and then extract the summed one-dimensional Lyα profile of a large (2.′′ 4 in diameter) region centered on this peak.We measure the wavelengths at which the Lyα emission from this large region meets the continuum on either side of the line, typically ∼ −900 to +1200 km s −1 , and use these as the wavelength limits of a 10. ′′ 5 × 10. ′′ 5 (i.e.35 × 35 0. ′′ 3 pixels) subcube centered on the Lyα 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 Lyα 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 Lyα subcube to create an emission-only cube.Finally, this cube is integrated along the wavelength axis to construct the continuum-subtracted Lyα surface brightness images shown in Figure 2. We note that, in general, accurate modeling and subtraction of the continuum underlying Lyα 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 Lyα 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 (Lyα) and black (continuum) contours in Figure 2.This comparison also shows that in most cases the Lyα and continuum peaks are spatially coincident.We measure the total Lyα flux of each extended halo by summing the largest connected region with S/N > 1 (corresponding to a surface brightness of ∼ 1-2 × 10 −18 erg s −1 cm −2 arcsec −2 for targets with approximately 5 hrs of integration), and calculate an effective circular radius r eff Lyα = (A/π) 1/2 , where A is the area of the region.The total Lyα fluxes range from 0.4 7×10 −16 erg s −1 cm −2 , corresponding to Lyα luminosities ranging from 1.6 × 10 42 to 2.6 × 10 43 erg s −1 with median 10 43 erg s −1 , while the radii vary between 16 and 30 kpc.
Lyα escape fractions are computed by comparing the total Lyα flux with the predicted flux calculated from the dustcorrected Hα emission4 and the theoretical Lyα/Hα ratio.Estimates of the Lyα escape fraction typically assume an intrinsic ratio (Lyα/Hα) int = 8.7 (e.g.Hayes 2015), but the Lyα/Hα 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 ∼ 100 -700 cm −3 , (Lyα/Hα) int ranges from 8.4 to 9.1. 5We therefore adopt an intrinsic ratio for each object that depends on the electron density; these ratios are listed in the column labeled (Lyα/Hα) 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 ∼ 450 cm −3 , corresponding to (Lyα/Hα) int = 8.9.
The resulting escape fractions range from 0.04 to 1.2, with a median of 0.22.We measure f Lyα esc = 1.23 ± 0.096 for Q2206-BX151, suggesting ∼ 100% of the Lyα emission escapes; however, as discussed in Section 2.1, the Hα 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 Lyα 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 Lyα esc = 0.96 ± 0.32 for Q0449-BX115, may also be overestimated due to the use of E(B − V ) cont rather than E(B − V ) neb for the dust correction; in practice, E(B − V ) cont = 0 for Q0449-BX115, so no dust correction was applied.
We next calculate the total Lyα equivalent width by dividing the total Lyα 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 ∼ 200 Å (for Q2206-BX151), roughly the maximum value expected from a normal stellar population (Charlot & Fall 1993).All Lyα measurements are given in Table 3.
Finally, we note that in this work we focus on the spectral properties of the extended Lyα emission in this strongly Lyα-emitting subset of the KBSS sample; an analysis of the structural and spectral properties of Lyα emission in the full KBSS-KCWI sample will be presented elsewhere.

THE SPECTRAL PROFILES OF SPATIALLY EXTENDED Lyα EMISSION
In this section we analyze the spatial variations of the spectral Lyα 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 twodimensional 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 Lyα 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 ≥ 3.In practice, most (∼ 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 ≥ 3 shown and each bin indicated by a different color.
Once the bins are identified, we extract the spectrum and variance of each bin from the continuum-subtracted Lyα datacube, using the mean of the individual spaxels for bins con-sisting of more than one spaxel.We then fit the Lyα profile of each bin with a double asymmetric Gaussian function in velocity space, defined as where A blue and A red and v 0,blue and v 0,red are the amplitudes and peak velocities of the blue and red components respectively.The asymmetric line width σ is further defined as σ = 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 Lyα 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 Lyα 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 Lyα peak separation for each spaxel as ∆v peak = v 0,red − v 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 timeconsuming 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 Lyα flux and halo size.The asymmetric Gaussian fits to each spaxel are overplotted in orange.
Maps of the Lyα 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 ∼ 0.1-1 and peak separations ranging from ∼ 300 to ∼ 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.
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 Lyα surface brightness peak.This colorcoding 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σ, 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σ 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, Lyα 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σ) 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σ 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σ 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 Lyα peak properties, Q0142-BX186, is the faintest object in the sample, with total Lyα 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.

SPATIALLY AVERAGED Lyα PROFILES
In order to determine general trends and study the spectral properties of the Lyα 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.

Annular Lyα profiles
We first study the average variation of the Lyα 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 Lyα images, we bin each halo in single spaxel (0. ′′ 3) radial increments.The spectra of all the spaxels in each bin are then summed, and the resulting Lyα 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 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 ×, the red peak is undetected for spaxels marked with a +, and the color of these spaxels indicates the 1σ upper limit in the case of blue non-detections and the 1σ 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 Lyα and continuum surface brightnesses respectively, as in Figure 2. ) Figure 6.Maps of the Lyα 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 ×, the red peak is undetected for spaxels marked with a +, and the color of these spaxels indicates the 1σ upper limit in the case of blue non-detections and the 1σ 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 Lyα and continuum surface brightnesses respectively, as in Figure 2.  measuring f tr , defined as the fraction of the total emission within ±100 km s −1 of the trough. 6 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 ∼ 0.2 ± 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 ∼ 16 kpc have flux ratios > 0.6 at that radius.The trough flux fraction f tr ranges from < 0 to ∼ 0.1 at the center of the halos, and rises consistently with radius for most objects.The largest measured values are f tr ∼ 0.2, found in the outer halos of Q0207-BX87 and Q2343-BX418.

Q0142-BX165
Trends with the average peak separation are somewhat more complicated.Most (10/12) of the halos have a central peak separation of ∼ 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 ∼ −100 km s −1 such that peak separations in the outer halo are ∼ 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, 6 ftr differs slightly from the quantity fcen defined by Naidu et al. (2022), who measure the fraction of flux escaping within ±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 peak separation is closely related to the trough depth f tr , decreasing as f tr increases.

Gradients in Lyα peak ratio and separation
While the binned, annular profiles described above are useful to characterize general trends in the extended Lyα 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 Lyα 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 Lyα 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 Lyα 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 Lyα 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 Figure 10.Results of the line profile measurements of the annular spectra described in Section 5.1.Top row: Lyα peak ratio vs. radius; Lyα peak separation vs. radius; and peak separation vs. ratio.Middle row: Lyα peak ratio vs. normalized Lyα surface brightness; Lyα peak separation vs. surface brightness; and surface brightness vs. radius.Bottom row: The fraction of total flux within ±100 km s −1 of the trough between the peaks ftr vs. radius; ftr vs. peak separation.
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 Lyα 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 Lyα emission.

MODELING LYα EMISSION AND LOW IONIZATION INTERSTELLAR ABSORPTION LINES
In the previous sections we have shown that the spectral morphology of Lyα 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 Lyα 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.

Lyα radiative transfer modeling
To extract physical properties of the gas in the halos from the observed Lyα profiles, we perform Monte Carlo radiative transfer (MCRT) modeling of the Lyα line.In contrast to the majority of previous studies in which spatially integrated Lyα spectra are modeled, in this work we attempt to fully leverage the power of KCWI and reproduce the spatially varying trends of the observed Lyα profiles.
Following a similar methodology to Li et al. (2022), we model the spatially resolved Lyα profiles using the multiphase, clumpy model.Each model is a 3D spherically symmetric region that emulates a galactic halo with a Lyα emitting source located at its center and two phases of gas: cool (∼ 10 4 K) H I clumps and a hot (∼ 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 Lyα profiles, primarily by producing additional absorption near the Lyα line center. 7In reality, such a component may correspond to the low column density absorbers (log N HI ≲ 10 17 cm −2 ) that provide additional Lyα scatterings in a galactic outflow (see e.g.Section 7.3 of Dijkstra & Kramer 2012).After interacting with these two phases of gas, Lyα 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 Lyα profiles.
In practice, we construct a grid of multiphase, clumpy models for fitting the Lyα spectra by varying the five most important physical parameters: F V , the volume filling factor of the clumps; log N HI, cl , the H I column density of the clumps; σ cl , the random velocity dispersion of the clumps; v cl , the radial outflow velocity of the clumps; and log n HI, ICM , the residual H I number density in the ICM. 8An additional parameter, ∆v, is used in post-processing to determine the deviation between the best-fit systemic redshift of the Lyα 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 V corresponds to a cloud covering factor f Lyα c (the average number of clumps per line of sight) of ∼ 1 − 10, which is similar to or moderately larger than the critical threshold f Lyα c, crit .Here f Lyα c, crit 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 Lyα c ≃ f Lyα c, crit ) where the Lyα RT in a multiphase, clumpy medium does not fully converge to the homogeneous shell model (Gronke et al. 2016(Gronke et al. , 2017;;Li & Gronke 2022).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 Lyα halos.Specifically, the momentum equation of an H I clump can be written as (Murray et al. 2005;Martin 2005): 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 −α .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 α less than 2. Assuming the gravitational potential is an isothermal sphere, we have M(r) = 2σ 2 cl r/G, where σ cl is the velocity dispersion of the clumps.Equation 2 can then be analytically solved as: where r min is the inner cutoff (or "launching") radius that satisfies v(r min ) = 0, and v cl, ∞ = 2Ar 1−α min /(α − 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 cl, ∞ due to the gravitational deceleration term; even the maximum radial v(r) is usually several hundred km s −1 smaller than v cl, ∞ .Following Dijkstra & Kramer (2012), we have fixed α = 1.4 and left σ cl and v cl, ∞ as the free parameters in this model.We set r min , so that r rmin ∈ [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 r rmin , so the following analysis applies to Lyα halos of varying physical sizes.
For each multiphase, clumpy model on the grid, MCRT has been performed on 10 4 Lyα photon packages emitted at the center of the simulation sphere in the form of a normalized Gaussian intrinsic spectrum N (v, µ = 0, σ = σ i, cl ), where σ i, cl = 12.85 km s −1 is the canonical thermal velocity dispersion of the H I gas in the clumps at T = 10 4 K. 9 The H I clumps with a constant column density N HI, cl are placed uniformly radially, so that their number density n cl ∝ 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 Lyα profiles by separating all the photons into three spatial bins according to their last-scattering impact parameters: b/b max ∈ (0, 0.25], (0.25, 0.50] and (0.50, 0.75], where b max is the largest impact parameter of the scattered Lyα 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 max and the halo radius r h is negligible, and we simply We only include the photons within 75% of b max (or equivalently, within the inner ∼ 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 ≤ 0.33r max , 0.33r max < r ≤ 0.67r max , and 0.67r max < r ≤ r max , where r 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 max , so that b/b max ∈ (0, 1 3 ], ( 1 3 , 2 3 ] and ( 2 3 , 1].Our fitting pipeline employs the python nested sampling package dynesty (Skilling 2004(Skilling , 2006;;Speagle 2020).At each visited point of the parameter space, the pipeline executes the following three steps: (1) calculate three binned Lyα 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 interpolation10 of the flux densities of the adjacent grid model spectra at the corresponding wavelength.The three binned Lyα model spectra are then convolved with a Gaussian function with σ = 65 km s −1 (the KCWI line spread function [LSF]) to mimic the finite instrumental resolution; (2) compare each binned model spectrum to an observed Lyα 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σ confidence intervals) of the samples in the marginalized PDF.

Metal absorption line modeling
In addition to the Lyα 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 λ1260 and C II λ1334 transitions,11 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 Lyα emission, i.e.Equations 2 and 3 with two free parameters: the clump velocity dispersion σ cl and the asymptotic maximum clump outflow velocity v cl, ∞ .Assuming that the absorption lines are saturated12 (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 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: where r min is the launching radius and f c,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 Lyα 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 Lyα modeling in §6.1, we fix α = 1.4 in the clump radial velocity profile and set r min = 0.1 kpc with r rmin ∈ [1, 100].Note again that only the ratio r rmin (rather than r or r 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: σ cl , v cl, ∞ , f c, max and γ.We use flat priors for the fitted parameters (which can vary continuously): , and γ ∈ [0.1, 2.0].We restrict γ 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 ∼ 100 kpc, see Figure 21 of Steidel et al. 2010 andRudie et al. 2019).

Modeling results & interpretation
Our modeling of both the spatially resolved Lyα emission and the UV absorption lines has achieved the following principal results: (1) reproducing the radially varying, spatially resolved Lyα profiles; (2) reproducing the radial trends of several important physical quantities of the Lyα 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 Lyα emission and metal absorption lines.
We present the modeling results for the spatially resolved Lyα spectra and the average line profile of Si II λ1260 and C II λ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 Lyα spectra (black); the middle row and the first panel of the bottom row show a comparison between the radial trends of 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-σ uncertainties).Note that the impact parameters may be slightly different for the model and the data: the models are binned consistently as b/bmax ∈ (0, 1 3 ], ( 1 3 , 2 3 ] and ( 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-σ uncertainties shown in grey) of Si II λ1260 (blue) and C II λ1334 (orange) profiles, as well as a comparison of clump radial outflow velocity profiles inferred from Lyα 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.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 Lyα 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 Lyα 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 inte-grated vs. spatially resolved Lyα modeling, and we discuss caveats to the models in Section 6.3.4.

Radial trends
The modeling results show that our multiphase, clumpy model is able to reproduce the spatially resolved Lyα 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 Lyα 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 Lyα RT perspective, the peak separation, which reflects the most likely frequencies at which the Lyα photons escape, is directly related to the Lyα optical depth of the system.The optical depth, which is the product of the Lyα cross-section 13 and the H I column density of the absorber, can therefore be expressed as a function of the temperature 13 Strictly speaking, the peak separation is also related to the gas outflow velocity, since the Lyα 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 column density of the absorber. 14The blue-to-red peak flux ratio, however, is negatively correlated with the H I gas outflow velocity as seen by the Lyα 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 Lyα 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 14 For example, the peak separation of Lyα photons that escape from an opaque, static H I sphere is ∆v peak ≃ 320 (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 ′ 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.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;15 (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. 16These three effects are presumably responsible for the observed radial variation of the spatially binned Lyα 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 V , log N HI, cl , σ cl , v cl, ∞ , log n HI, ICM , ∆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 Lyα 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 cl (r) ≃ 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 Lyα 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.
Incidentally, our model has also reproduced the decreasing trend of Lyα 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 Model III : Higher nICM at outer halo Figure 14.Experiments designed to test our hypotheses for the differences between Lyα photons that escape at low and high impact parameters.In each of the four subpanels, three binned model Lyα spectra are shown according to their last-scattering impact parameters: b/rh ∈ (0, 1 3 ] (green solid), ( 13 , 2 3 ] (red dash-dotted) and ( 23 , 1] (blue dotted), where rh is the radius of the modeled halo.Left: The fiducial model with (FV, log NHI,cl, σcl, vcl,∞, log nHI,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 ncl(r) ≃ 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 Lyα spectral morphology (i.e.peak separation, peak flux ratio and trough flux fraction), hence supporting our explanation.
of the total halo area used for the spatially resolved Lyα modeling.This overall consistency adds further credence to our multiphase, clumpy RT model.

Best-fit parameters
One of the most interesting discoveries from our modeling is that the clump outflow velocities inferred from Lyα emission and the low ionization metal absorption lines can be mutually consistent, with typical values of ∼ 400 − 600 km s −1 obtained for both.The mismatch between the gas outflow velocities inferred from Lyα and from metal absorption lines has been a long-standing problem.For example, it is reported that the ≲ 150 km s −1 outflow velocities of the shell model required to match the Lyα profiles of local starburst and green pea galaxies are much lower than the ≳ 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 Lyα RT models has been little explored, possibly due to the belief that the Lyα 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 Lyα 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 Lyα profiles in the high outflow velocity regime (v cl, max ≳ 400 − 600 km s −1 ).
In our sample, consistency (accounting for uncertainties) between the clump outflow velocities inferred from Lyα and metal absorption lines is achieved in 8 / 12 objects.17Such a high success rate demonstrates the feasibility of matching both the observed Lyα 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 Lyα and absorptionline-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 Lyα results incorporate gas on the far side and at large impact parameters that is not seen in absorption.
The best-fit radial velocity profile of the clumps in the multiphase, clumpy model typically exhibits a rapid acceleration phase to v cl = v cl, max within 1 ≲ r rmin ≲ 10 followed by a gradual deceleration 18 (or v cl ≃ constant) phase at r rmin ≳ 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.(0.04, 17.5, −6.5, 125) and vcl,∞ = (500, 600, 700, 800, 900) are shown with different colors and linestyles.Upper: spatially integrated Lyα model spectra with different vcl,∞ values.As vcl,∞ 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 vcl,∞.Note that vcl(r) and vcl,∞ are positively correlated, but vcl(r) is always smaller than (typically by several hundred km s −1 ) vcl,∞ due to the effect of gravitational deceleration.
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, andreferences therein).It is possible, however, that the interaction between the hot phase and the Lyα photons is dominated by the decelerated, semi-static hot gas, as suggested by the deep troughs at line center in the observed Lyα profiles.A larger sample with more diverse Lyα 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 Lyα modeling, the best-fit clump volume filling factors (F V ) range from 0.06 to 0.16 (corresponding to ∼ 5 − 10 clumps on average per line-of-sight 19 ), and the best-fit clump column densities (N HI, cl ) range from ∼ 10 17.6 to 10 18.8 cm −2 .The total H I column densities (N HI, total ≃ 4 3 f cl N HI, cl = (r h /r cl ) F V N HI, cl , Gronke et al. 2016) of the bestfit models range from ∼ 10 18.5 to 10 19.9 cm −2 .Here N HI, total represents the inferred total H I column density of the modeled halo that a Lyα photon typically interacts with, either via scattering or free-streaming; the scattered, out-of-resonance Lyα 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 HI, ICM ≃ n HI, ICM r h ) range from ∼ 10 15 to 10 16 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 Lyα sources are mostly consistent with the systemic redshifts determined from nebular emission lines (|∆v| < 50 km s −1 for 8 / 12 objects).The best-fit clump velocity dispersions (σ cl ) are all smaller than 150 km s −1 and span a similar range to the observed nebular emission line widths (∼ 50 − 120 km s −1 ).We compared the best-fit σ cl values with the MOSFIRE Hband ([O III] and Hβ) and K-band (Hα) 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 dispersions 20 are all smaller than 75 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 ≳ 500 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 c,max range from ∼ 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 (γ) range from ∼ 1.1 to 1.7, corresponding to a mass-conserving-like (or more gradual) decrease in the number density of the clumps.A γ 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, γ will be exactly 2 due to the geometric volume increase at large radii (dV ∝ 4π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 cl, max inferred from the Lyα 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.

Advantages of spatially resolved Lyα modeling
In this section, we demonstrate the advantages of spatially resolved Lyα modeling by comparing it to the spatially integrated Lyα modeling that has typically been carried out in previous works.Assuming that for a Lyα-emitting source of interest, only a spatially integrated Lyα 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 Lyα 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 Lyα 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 Lyα emission.We model scenario (2) with spatially integrated models that include the photons with b/r h ≲ 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 σ cl and v cl, ∞ 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 Lyα 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 Lyα 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 Lyα 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 HI and n 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 HI from spatially integrated modeling of the central region are larger on average by a factor of 1.5, while n 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 HI, ICM higher than the maximum value allowed by the model grid.The overestimation of n 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 Lyα 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 Lyα photons, or fail to account for the radial variation of the Lyα morphological properties.In comparison, spatially resolved Lyα modeling has the advantage of fully leveraging the spatial variation in the Lyα 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.

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 Lyα photons will eventually escape from the simulation region and contribute to the emergent model spectra.Considering that the actual Lyα 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 Lyα photons is representative of the Lyα 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 Lyα 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 Lyα 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 10 6 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.

SUMMARY AND DISCUSSION
We have presented KCWI integral field spectroscopy and radiative transfer modeling of spatially extended Lyα emission in a sample of 12 relatively low mass (M ⋆ ∼ 10 9 M ⊙ ), 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 Lyαemitters.The sample galaxies have specific star formation rates ∼ 4 times larger than than that of their z ∼ 2 parent sample, and may more closely resemble galaxies at earlier epochs of cosmic history.Our primary results are as follows: 1.All of the galaxies show strong, double-peaked Lyα emission (see Section 3 and Figure 1) and spatially extended Lyα halos, with luminosities ranging from 3 × 10 42 to 3 × 10 43 erg s −1 and radii between 16 and 30 kpc (Figure 2).
2. We fit double asymmetric Gaussian profiles to the Lyα 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 ∆v 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-tored 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. We also construct spatially averaged Lyα 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 tr , the fraction of total flux escaping within ±100 km s −1 of the trough between the peaks (Figure 10).The blueto-red flux ratio increases consistently with radius for most objects in the sample, with a typical central value of ∼ 0.2; all objects that can be measured at a radius ≳ 16 kpc have peak flux ratios > 0.6 at that radius.f tr also increases with radius for most of the sample.
Trends with peak separation are more complex, but the typical central peak separation is ∼ 600 km s −1 , with a moderate decrease toward the outer halo.
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 Lyα blue-to-red flux ratio, peak separation, trough flux fraction ftr 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: Lyα blue-to-red flux ratio, peak separation, and ftr 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.
halos have sightlines for which the peak ratio increases (typically from ∼ 0.2 to ∼ 1) or the peak separation decreases (typically from ∼ 600 -700 to ∼ 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. Using a new suite of Lyα radiative transfer simulations, we model the spatially resolved Lyα 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 .The clumps reach a typical maximum velocity of ∼ 500 km s −1 and have H I column densities of ∼ 10 17.6 to 10 18.8 cm −2 , while the total N HI of the best-fit models ranges from ∼ 10 18.5 to 10 19.9 cm −2 .The clumps are embedded in a hot interclump medium with residual N HI, ICM ∼ 10 15 -10 16 cm −2 .Best-fit parameters of the models are given in Table 5.
6.We find that the trend in Lyα 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 be-fore they escape due to the decrease in clump covering fraction with radius.The Lyα 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 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. We self-consistently model the mean low ionization absorption line profile of each object, employing the same radially varying velocity model used for the Lyα 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 ≳ 500 km s −1 , in broad agreement with the velocities inferred from Lyα; exact matches may not be expected because the down-thebarrel UV spectra and the radially binned Lyα emission are not probing entirely the same regions of the halos.This agreement alleviates a long-standing discrepancy between outflow velocities inferred from Lyα shell models and the UV absorption lines.
8. Finally, we compare the results of the spatially resolved Lyα modeling with those obtained from applying the same model to single, spatially integrated Lyα 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 Lyα 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 ∼ 2 starforming galaxies: a multiphase, clumpy medium in which cool (∼ 10 4 K), outflowing gas clumps are embedded in a hot (∼ 10 6 K), highly ionized, diffuse medium with low-density residual H I. The clumps typically have H I column densities of ∼ 10 18 cm −2 and provide a total column density of ∼ 10 19 cm −2 , and the Lyα 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 ∼ 100 km s −1 , and are accelerated to high radial outflow velocities of ≳ 500 km s −1 at large impact parameters, which give rise to both the asymmetric Lyα profiles and broad low ionization metal absorption lines.The hot ICM is nearly static and has a low total H I column density (∼ 10 15 -10 16 cm −2 ), but is essential to shaping the emergent double-peaked Lyα profiles as it provides additional scattering that produces the absorption trough at line center.

Central Lyα profiles and LyC escape
With this physical model of the CGM in mind, we revisit the spatially integrated central Lyα profiles shown in Figure 1 and assess how (or if) quantities measured from these profiles relate to the properties of the CGM inferred from the spatially resolved modeling; such a comparison may aid in the interpretation of Lyα profiles when information from the outer halo is unavailable.In Figure 16 we compare the total N HI and n HI, ICM from the models with the peak ratio, peak separation and trough flux fraction f tr and the mean low ionization absorption equivalent width W LIS measured from the spatially integrated one-dimensional spectra, as well as the equivalent width vs. the Lyα profile properties in the bottom row.Darker points indicate the results of the spatially resolved Lyα modeling, while the fainter points are the result of modeling the central spatially integrated profiles (scenario 1 in Section 6.3.3).The lower corner of each panel gives the p-value resulting from a Spearman correlation test, with values of p < 0.1 highlighted in red.While none of the correlations are formally (> 3σ) significant, the strongest trends (∼ 2.75-3σ) relate to the H I density in the ICM, which tends to be higher for larger peak separations, lower f tr , and larger low ionization equivalent width.We also find that smaller peak separations and higher values of f tr tend to be associated with lower W LIS .All correlations involving the total N HI or the blue-to-red peak ratio have significance levels ≤ 1.3σ.
These results broadly support our conclusion in Section 6.3.1 that the trough flux fraction can be understood as an indication of low N HI in the ICM.Note, however, that the potential relationship between the central f tr and modeled n HI, ICM relies on the results inferred from spatially resolved modeling of the extended halo; modeling the central profiles alone results in significantly higher values of n HI, ICM , half of which are higher than the upper limit of the current model grid.
Previous work has suggested that significant Lyα flux at the systemic velocity may be an indication of LyC escape (e.g.Rivera-Thorsen et al. 2019;Naidu et al. 2022); if ionizing photons emerge through optically thin channels between clumps, then the transparency of the ICM is a key property governing LyC escape.A low covering fraction of neutral gas and significant residual intensity in the low ionization absorption lines are also likely related to LyC escape (e.g.Heckman et al. 2011;Reddy et al. 2016;Chisholm et al. 2018), so the potential relationship between n HI, ICM and W LIS is also unsurprising.
Given the results of the spatially resolved Lyα modeling, we expect the peak separation to be most closely related to the total H I column density; however, there is no significant correlation between the central peak separation and N HI, total from the spatially resolved models.This lack of correlation may be due to the small sample size and the lack of dynamic range in peak separation, as 10 of the 12 objects in the sample have central peak separations between 500 and 700 km s −1 .These peak separations are also larger than the ∼ 200-500 km s −1 which ∆v peak is observed to correlate with the LyC escape fraction in local galaxies (Izotov et al. 2021).We do observe potential relationships between the peak separation and both n HI, ICM and W LIS ; these may be due to the strong correlation between the peak separation and f tr .Modeling of a larger sample with a wider range of central peak separations will clarify the relationship between ∆v peak and N HI, total .
There are no observations covering wavelengths below the Lyman break for the galaxies in our sample, so we have no constraints on their LyC emission.However, based on the criteria discussed above involving the peak separation or central flux fraction, we would not expect most of the galaxies in the sample to have significant LyC emission.Possible exceptions are the two most likely LyC candidates, Q0821-MD36 and Q0207-BX87, which have the highest trough flux fractions f tr ∼ 0.1, relatively narrow peak separations, and the second and third highest Lyα equivalent widths in the sample (after Q2206-BX151).

Future prospects
Although the inclusion of spatially resolved information increases the power of the radiative transfer modeling, we are still limited by the assumption of symmetry: we fit radially binned spectra with spherically symmetric models, but as we have shown, real halos show significant azimuthal variation (Figure 11).However, insights obtained from the modeling can aid in the interpretation of the variations across a given halo, at least qualitatively.Because the increase in blue-tored peak ratio with radius is largely a geometric effect due to the decrease in the line-of-sight component of the outflow velocity, portions of the halos for which there is little change in the peak ratio with radius likely correspond to regions for which the velocity still has a significant component along the observer's line of sight even in the outskirts of the halo.More broadly, azimuthal variations in the peak ratio are indicative of velocity asymmetries and non-radial gas motions at large radii.Similarly, variations in the peak separation in the outer halo suggest varying H I column densities in the CGM, with regions for which ∆v peak does not decrease with radius likely having higher N HI .Future modeling that does not assume azimuthal symmetry is needed in order to quantify these conclusions.
Finally, while the objects in this sample are likely to be more typical of galaxies at higher redshifts than of the general z ∼ 2 population, extending the analysis of double-peaked Lyα profiles to more distant galaxies will be challenging.For example, the median redshift of the MUSE sample studied by Leclercq et al. (2020) is z = 3.8, while that of our KCWI sample is z = 2.3, and this difference in redshift results in a median decrease in surface brightness of a factor of 4.5 for the higher redshift sample.In addition, the blue-to-red Lyα peak ratio decreases with increasing redshift due to Lyα absorption by the IGM (Laursen et al. 2011;Hayes et al. 2021), and the mean IGM transmission of Lyα drops strongly from ≳ 80% at z ≈ 2.3 to ∼ 45% at z ≈ 3.8 (Rudie et al. 2013;Inoue et al. 2014).The combination of these effects results in a typical factor of ≳ 6 decrease in the surface brightness of the blue peak at z = 3.8 relative to z = 2.3.These effects will of course be even more significant at z > 4.
Given the power of the double-peaked Lyα profile for constraining the kinematics and column density of the CGM, we therefore expect that integral field observations of galaxies at z ∼ 2 will only grow in importance.As new observations from the James Webb Space Telescope precisely measure the properties of galaxies at both z ∼ 2 and in the reionization era, it will be increasingly possible to robustly identify z ∼ 2 analogs of reionization-era sources and quantify their CGM via spatially extended Lyα emission.
The authors thank the referee for a detailed and constructive report, and Crystal Martin, Peng Oh and Tim Heckman for helpful discussions.D.K.E. is supported by the US National Science Foundation (NSF) through Astronomy & Astrophysics grant AST-1909198. C.C.S, Z.L. and Y.C. are supported in part by NSF grant AST-2009278.R.F.T. is a Cottrell Scholar receiving support from the Research Corporation for Science Advancement under grant ID 28289 and from the Pittsburgh Foundation under grant ID UN2021-121482.Some data presented herein were obtained at the W. M. Keck Observatory from telescope time allocated to the National Aeronautics and Space Administration through the agency's scientific partnership with the California Institute of Technology and the University of California.The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.The authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Maunakea has always had within the indigenous Hawaiian community.We are most fortunate to have the opportunity to conduct observations from this mountain.D.K.E.acknowledges in Milwaukee that we are on traditional Potawatomi, Ho-Chunk and Menominee homeland along the southwest shores of Michigami, North America's largest system of freshwater lakes, where the Milwaukee, Menominee and Kinnickinnic rivers meet and the people of Wisconsin's sovereign Anishinaabe, Ho-Chunk, Menominee, Oneida and Mohican nations remain present.

APPENDIX A. MODELING RESULTS OF THE FULL SAMPLE
In Figures 17 to 22 below we present the results of Lyα and metal absorption line modeling for all objects except Q0207-BX144 (already shown in Figure 12).In each panel, the top row shows the best-fit RT models (red) to the spatially resolved Lyα 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 Lyα emission and the average metal absorption line.Software: KCWI DRP 21 , astropy (Astropy Collaboration et al. 2013, 2018), dynesty (Skilling 2004(Skilling , 2006;;Speagle 2020), spectral-cube (Ginsburg et al. 2019), SExtractor (Bertin & Arnouts 1996), cmasher (van der Velden 2020), seaborn (Waskom et al. 2020)

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

Figure 2 .
Figure 2. Continuum-subtracted Lyα images.White contours show the Lyα surface brightness, with the same levels in each panel: 1 × −18 (dotted), 5 × 10 −18 (dashed), 1 × 10 −17 (dash-dot), and 2 × 10 −17 (solid) erg s −1 cm −2 arcsec −2 , and black contours indicate the adjacent UV continuum measured in a rest-frame 75 Å window redward of the Lyα emission line.ples remain incomplete at the low masses characteristic of our targets).The [O III]/[O II] ratios 1 provide an estimate of the degree of excitation in the H II regions.As expected given the sample selection criterion of high [O III]/Hβ ratios, our targets fall at the upper end of the O32 distribution for the

Figure 3 .
Figure 3. Representative results of the Voronoi binning of the Lyα surface brightness images.The bins include only individual pixels with S/N ≥ 1, and are constructed to achieve a target S/N ≥ 3. Bins generally consist of single pixels, except in the outer halo where each object has a few bins of ∼ 2-5 pixels.Contours showing Lyα surface brightness are the same as in Figure 2, and show that regions reaching the target S/N generally have surface brightness ≳ 5 × 10 −18 erg s −1 cm −2 arcsec −2 (dashed contour).

Figure 4 .
Figure 4.The Lyα 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 Lyα surface brightness of the full modeled region of Q0449-BX110 on a logarithmic scale, with the portion shown here outlined in orange.The Lyα profiles are shown in the order of their corresponding spaxels in the orange box.Contours of Lyα surface brightness are the same as in Figure 2.

Figure 5 .
Figure5.Maps of the Lyα 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 ×, the red peak is undetected for spaxels marked with a +, and the color of these spaxels indicates the 1σ upper limit in the case of blue non-detections and the 1σ 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 Lyα and continuum surface brightnesses respectively, as in Figure2.

Figure 7 .
Figure 7. Lyα 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 Lyα surface brightness image.Spearman r and p values are given in the upper right of each panel, with significant (> 3σ) correlations labeled in red.With p = 2.83 × 10 −3 , Q1700-BX729 has 2.99σ significance.

Figure 8 .
Figure 8. Significance levels from pairwise Spearman correlation tests of distance from the peak of the Lyα surface brightness image, Lyα surface brightness, peak ratio and peak separation.

Figure 9 .
Figure9.Normalized annular Lyα profiles, constructed by binning all spaxels that have S/N > 2 in the continuum-subtracted Lyα images in single spaxel (0. ′′ 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.

Figure 11 .
Figure11.Top two rows: Blue-to-red Lyα 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 Lyα peak separation, with maximum gradients indicated by gold squares and minimum gradients by light yellow diamonds.The annular averages from Figure10are also shown as red stars in all panels.See Section 5.2 for details.

Figure 12 .
Figure12.Modeling results of the annular-averaged, spatially resolved Lyα spectra and of the average line profile of Si II λ1260 and C II λ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 Lyα spectra (black, with 1-σ 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-σ uncertainties).Note that the impact parameters may be slightly different for the model and the data: the models are binned consistently as b/bmax ∈ (0, 1 3 ], ( 1 3 , 2 3 ] and ( 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-σ uncertainties shown in grey) of Si II λ1260 (blue) and C II λ1334 (orange) profiles, as well as a comparison of clump radial outflow velocity profiles inferred from Lyα 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.

Figure 13 .
Figure13.Schematic of the escape of two Lyα 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 Lyα 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 ′ 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 ′ 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 cl, ∥ < v ′ cl, ∥ , 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 ′ 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.

FFigure 15 .
Figure 15.Examples of Lyα model spectra with different clump outflow velocities showing the pattern in the change of the blue-tored peak flux ratio.Five models with (FV, log NHI,cl, log nICM, σcl) =(0.04, 17.5, −6.5, 125)  and vcl,∞ = (500, 600, 700, 800, 900) are shown with different colors and linestyles.Upper: spatially integrated Lyα model spectra with different vcl,∞ values.As vcl,∞ 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 vcl,∞.Note that vcl(r) and vcl,∞ are positively correlated, but vcl(r) is always smaller than (typically by several hundred km s −1 ) vcl,∞ due to the effect of gravitational deceleration.

Table 2 .
Nebular Emission Line Measurements

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