The KBSS–KCWI survey: the connection between extended Ly α haloes and galaxy azimuthal angle at z ∼ 2–3 Free

ABSTRACT

We present the first statistical analysis of kinematically resolved, spatially extended |$\rm Ly\alpha$| emission around z = 2–3 galaxies in the Keck Baryonic Structure Survey (KBSS) using the Keck Cosmic Web Imager (KCWI). Our sample of 59 star-forming galaxies (z_med = 2.29) comprises the subset with typical KCWI integration times of ∼5 h and with existing imaging data from the Hubble Space Telescope and/or adaptive optics-assisted integral field spectroscopy. The high-resolution images were used to evaluate the azimuthal dependence of the diffuse |$\rm Ly\alpha$| emission with respect to the stellar continuum within projected galactocentric distances of ≲30 proper kpc. We introduce cylindrically projected 2D spectra (CP2D) that map the averaged |$\rm Ly\alpha$| spectral profile over a specified range of azimuthal angle, as a function of impact parameter around galaxies. The averaged CP2D spectrum of all galaxies shows clear signatures of |$\rm Ly\alpha$| resonant scattering by outflowing gas. We stacked the CP2D spectra of individual galaxies over ranges of azimuthal angle with respect to their major axes. The extended |$\rm Ly\alpha$| emission along the galaxy principal axes is statistically indistinguishable, with residual asymmetry of ≤2 per cent (∼2σ) of the integrated |$\rm Ly\alpha$| emission. The symmetry implies that the |$\rm Ly\alpha$| scattering medium is dominated by outflows in all directions within 30 kpc. Meanwhile, we find that the blueshifted component of |$\rm Ly\alpha$| emission is marginally stronger along galaxy minor axes for galaxies with relatively weak |$\rm Ly\alpha$| emission. We speculate that this weak directional dependence of |$\rm Ly\alpha$| emission becomes discernible only when the |$\rm Ly\alpha$| escape fraction is low. These discoveries highlight the need for similar analyses in simulations with |$\rm Ly\alpha$| radiative transfer modelling.

1 INTRODUCTION

Extended ‘haloes’ of diffuse Lyman-alpha (Ly α) emission, extending to many times larger radii than starlight, are nearly ubiquitous around rapidly star-forming galaxies at redshift z > 2 (e.g. Steidel et al. 2011; Wisotzki et al. 2016; Leclercq et al. 2017; Wisotzki et al. 2018). While there is not yet consensus on the dominant physical mechanism giving rise to diffuse |$\rm Ly\alpha$| haloes, there is general agreement on the list of potential sources, all of which depend on substantial cool hydrogen gas in the circumgalactic medium (CGM; e.g. Ouchi, Ono & Shibuya 2020): (1) resonant scattering of |$\rm Ly\alpha$| produced by recombination of gas photoionized by massive stars or active galactic nuclei (AGNs), i.e. a central source, where |$\rm Ly\alpha$| photons are subsequently scattered until they find optically thin channels to escape; (2) in situ photoionization of H i by the metagalactic ultraviolet (UV) ionizing radiation field combined with local sources and followed by recombination (sometimes called ‘fluorescence’); (3) accreting gas losing energy via collisional excitation of |$\rm Ly\alpha$| (sometimes referred to as ‘gravitational cooling’); and (4) emission from unresolved satellite galaxies in the haloes of larger central galaxies. In principle, the observed surface brightness (SB), spatial distribution, and kinematics of |$\rm Ly\alpha$| emission can discriminate between the various mechanisms and, perhaps more importantly, can provide direct information on the degree to which gas in the CGM is accreting, outflowing, or quiescent. |$\rm Ly\alpha$| emission from the CGM, if interpreted correctly, can provide a detailed map of the cool component of the dominant baryon reservoir associated with forming galaxies, as well as constraints on large-scale gas flows that are an essential part of the current galaxy formation paradigm.

Because a |$\rm Ly\alpha$| photon is produced by nearly every photoionization of hydrogen, the intrinsic |$\rm Ly\alpha$| luminosity of a rapidly star-forming galaxy can be very high and thus easily detected (Partridge & Peebles 1967). However, due to its very high transition probability, |$\rm Ly\alpha$| is resonantly scattered until the last scattering event gives it an emitted frequency and direction such that the optical depth remains low along a trajectory that allows it to escape from the host galaxy. When the |$\rm Ly\alpha$| optical depth is high in all directions, the vastly increased effective path length – due to large numbers of scattering events during the time the photon is radiatively trapped – increases the probability that the photon is destroyed by dust or emitted via two-photon mechanism. But |$\rm Ly\alpha$| photons that are not absorbed by dust grains or converted to two-photon radiation must eventually escape, with the final scattering resulting in the photon having a frequency and direction such that the photon can freely stream without further interaction with a hydrogen atom. The radiative transfer of |$\rm Ly\alpha$| thus depends in a complex way on the distribution, clumpiness, and kinematics of H i within the host galaxy, as well as on where and how the photon was produced initially. But the added complexity is counter-balanced by the availability of a great deal of information about the scattering medium itself that is otherwise difficult or impossible to observe directly: i.e. the neutral hydrogen distribution and kinematics in the CGM.

Since the commissioning of sensitive integral-field spectrometers on large ground-based telescopes – the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) on the Very Large Telescopes (VLT) of the European Southern Observatory and the Keck Cosmic Web Imager (KCWI; Morrissey et al. 2018) at the Keck Observatory – it has become possible to routinely detect diffuse |$\rm Ly\alpha$| emission haloes around individual galaxies at high redshift (e.g. Wisotzki et al. 2016; Leclercq et al. 2017) and to simultaneously measure the spatially resolved |$\rm Ly\alpha$| kinematics (e.g. Erb, Steidel & Chen 2018; Claeyssens et al. 2019; Leclercq et al. 2020). Such observations can then be interpreted in terms of simple expectations based on |$\rm Ly\alpha$| radiative transfer; for example, a generic expectation is that most |$\rm Ly\alpha$| emission involving scattered photons (i.e. those that must pass through an H i gas distribution before escaping their host) will exhibit a ‘double-peaked’ spectral profile, where the relative strength of the blueshifted and redshifted peaks may be modulated by the net velocity field of the emitting gas. In the idealized case of a spherical shell of outflowing (infalling) gas, one predicts that an external observer will measure a dominant red (blue) peak (Verhamme, Schaerer & Maselli 2006). The fact that most (∼90 per cent) of star-forming galaxy ‘down the barrel’ (DTB) spectra with |$\rm Ly\alpha$| in emission in the central portions exhibit dominant red peaks (e.g. Pettini et al. 2001; Steidel et al. 2010; Kulas et al. 2012; Trainor et al. 2015; Verhamme et al. 2018; Matthee et al. 2021) has led to the conclusion that outflowing gas dominates |$\rm Ly\alpha$| radiative transfer, at least at small galactocentric distances.

Essentially every simulation of galaxy formation (Faucher-Giguère et al. 2010) predicts that gaseous accretion is also important – particularly at high redshifts (z ≳ 2) – and this has focused attention on systems in which a double |$\rm Ly\alpha$| profile with a blue-dominant peak is observed, often cited as evidence for ongoing accretion of cool gas (e.g. Martin et al. 2014, 2016; Vanzella et al. 2017; Ao et al. 2020). Quantitative predictions of |$\rm Ly\alpha$| emission from accreting baryons depend sensitively on the thermal state and the small-scale structure of the gas (e.g. Faucher-Giguère et al. 2010; Goerdt et al. 2010; Kollmeier et al. 2010), leading to large uncertainties in the predictions. The role played by ‘local’ sources of ionizing photons over and above that of the metagalactic ionizing radiation field is likely to be substantial for regions near QSOs (Cantalupo et al. 2014; Borisova et al. 2016; Cai et al. 2019; O’Sullivan et al. 2020) but much more uncertain for star-forming galaxies, where the escape of scattered |$\rm Ly\alpha$| photons is much more likely than that of ionizing photons.

Models of |$\rm Ly\alpha$| radiative transfer have attempted to understand the dominant physics responsible for producing |$\rm Ly\alpha$| haloes around galaxies. Using photon-tracing algorithms with Monte Carlo simulations, Verhamme et al. (2006), Dijkstra (2014), Gronke & Dijkstra (2016), and Gronke et al. (2016) have explored the effects of resonant scattering on the emergent central |$\rm Ly\alpha$| line profile using various idealized H i geometries and velocity fields; in most cases, simple models can be made to fit the observed 1D profiles (Gronke 2017; Song, Seon & Hwang 2020). There have also been attempts to model or predict spatially resolved |$\rm Ly\alpha$|⁠, which almost certainly depends on a galaxy’s immediate environment (Zheng et al. 2011; Kakiichi & Dijkstra 2018), including both outflows and accretion flows, as well as the radiative transfer of |$\rm Ly\alpha$| photons from the site of initial production to escape (e.g. Faucher-Giguère et al. 2010; Lake et al. 2015; Smith et al. 2019; Byrohl et al. 2020). However, the conclusions reached as to the dominant process responsible for the extended |$\rm Ly\alpha$| emission have not converged, indicating that more realistic, high-resolution, cosmological zoom-in simulations may be required to capture all of the physical processes. Despite the variety of |$\rm Ly\alpha$| radiative transfer models to date, as far as we are aware, no specific effort has been made to statistically compare full 2D model predictions (simultaneous spatial and kinematic) to observed |$\rm Ly\alpha$| haloes.

Some insight into the relationship between galaxy properties and the kinematics and spatial distribution of cool gas in the CGM has been provided by studies at lower redshifts, where galaxy morphology is more easily measured. Statistical studies using absorption line probes have clearly shown that the strength of low-ionization metal lines such as Mg ii, Fe ii depends on where the line of sight passes through the galaxy CGM relative to the projected major axis of the galaxy – the ‘azimuthal angle’ (e.g. Bordoloi et al. 2011; Bouché et al. 2012; Kacprzak, Churchill & Nielsen 2012; Nielsen et al. 2015; Lan & Mo 2018; Martin et al. 2019), and the inclination of the galaxy disc relative to the line of sight (e.g. Steidel et al. 2002; Kacprzak et al. 2011). More recently, clear trends have also been observed for high ions (O vi) (Kacprzak et al. 2015). In general, these trends support a picture of star-forming galaxies in which high-velocity, collimated outflows perpendicular to the disc are responsible for the strongest absorption lines in both low and high ions, with low ions also being strong near the disc plane. Theoretically at least, accretion flows might also be quasi-collimated in the form of cold streams of gas that would tend to deposit cool gas near the disc plane (see e.g. Tumlinson, Peeples & Werk 2017). It is less clear how such a geometry for gas flows in the CGM would manifest as emission in a resonantly scattered line like |$\rm Ly\alpha$|⁠. One might expect that |$\rm Ly\alpha$| photons would escape most readily along the minor axis, since the large velocity gradients and lower H i optical depths of outflowing material both favour |$\rm Ly\alpha$| escape from the host galaxy. This picture is consistent with Verhamme et al. (2012), who showed that |$\rm Ly\alpha$| escape is enhanced when the simulated galaxies are viewed face-on.

Observations of low-redshift, spatially resolved |$\rm Ly\alpha$| emission have so far been limited to small samples – e.g. in the local universe (z < 0.2), using the Hubble Space Telescope (HST), the ‘|$\rm Ly\alpha$| reference sample’ (LARS; Östlin et al. 2014) has obtained images probing |$\rm Ly\alpha$| emission around galaxies in great spatial detail, reaffirming the complex nature of |$\rm Ly\alpha$| radiative transfer and its relation to the host galaxies. In most cases, LARS found evidence that extended |$\rm Ly\alpha$| emission is most easily explained by photons produced by active star formation that then diffuse into the CGM before a last scattering event allows escape in the observer’s direction. Although there are small-scale enhancements associated with outflows, even for galaxies observed edge-on the |$\rm Ly\alpha$| emission is perhaps smoother than expected (Duval et al. 2016).

At z > 2, |$\rm Ly\alpha$| emission is more readily observed but detailed analyses are challenged by the relatively small galaxy sizes (both physical and angular) and the need for high spatial resolution to determine the morphology of the stellar light. In this paper, we present a statistical sample of z > 2 galaxies drawn from a survey using KCWI of selected regions within the Keck Baryonic Structure Survey (KBSS; Rudie et al. 2012; Steidel et al. 2014; Strom et al. 2017). Since the commissioning of KCWI in late 2017, we have obtained deep integral field unit (IFU) data (∼5-h integrations) for >100 KBSS galaxies with z = 2–3.5, so that the |$\rm Ly\alpha$| line is covered within the KCWI wavelength range; some initial results from the survey have been presented by Erb et al. (2018) and Law et al. (2018). The 59 galaxies included in our current analysis are those that, in addition to the KCWI data, have also been observed at high spatial resolution by either HST or adaptive optics-assisted near-IR spectroscopy using Keck/OSIRIS. The overarching goal of the study is to evaluate the spatial and spectral distribution of |$\rm Ly\alpha$| emission within ≃5 arcsec as compared to the principle axes defined by the galaxy morphology on smaller angular scales by each galaxy’s UV/optical continuum emission. In particular, we seek to use the observed kinematics and spatial distribution of |$\rm Ly\alpha$| emission to evaluate whether the cool gas in the CGM of forming galaxies shows evidence for directional dependence – e.g., inflows or outflows along preferred directions – with respect to the central galaxy.

This paper is organized as follows. In Section 2, we describe the KBSS–KCWI sample; Section 3 introduces the high-resolution imaging and IFU data set; Section 4 covers the details on the measurement of the galaxy principle axes providing the definition of the galactic azimuthal angle; Section 5 presents results on the connection between Ly α haloes and galactic azimuthal angle; Section 6 looks into the connection between the |$\rm Ly\alpha$| azimuthal asymmetry and the overall |$\rm Ly\alpha$| emission properties; Section 7 checks higher-order azimuthal asymmetry of |$\rm Ly\alpha$| emission by dividing the sample into finer azimuthal bins; and Section 8 discusses the implications of the results, with a summary in Section 9. Throughout the paper, we assume a ΛCDM cosmology with Ω_m = 0.3, |$\Omega _\Lambda =0.7$|⁠, and h = 0.7. Distances are given in proper units, i.e. physical kpc (pkpc).

2 THE KBSS–KCWI GALAXY SAMPLE

In late 2017, we began using the recently commissioned KCWI (Morrissey et al. 2018) on the Keck II 10-m telescope to target selected regions within the survey fields of the KBSS (Rudie et al. 2012; Steidel et al. 2014; Strom et al. 2017). The main goal of the KCWI observations has been to detect diffuse emission from the CGM (within impact parameter, D_tran ≲ 100 pkpc) of a substantial sample of rapidly star-forming galaxies and optically faint AGN host galaxies, reaching SB sensitivity of ∼5 × 10⁻²⁰ ergs s⁻¹ cm⁻² arcsec⁻² (1σ) for unresolved emission lines. Such limiting SB would allow detection of the extended |$\rm Ly\alpha$| haloes of individual galaxies at redshifts z ∼ 2–3 (e.g. Steidel et al. 2011) and would be capable of detecting extended diffuse UV metallic cooling line emission as predicted by simulations of galaxies with comparable mass and redshift (e.g. Sravan et al. 2016).

KCWI offers three selectable image slicer scales, each of which can be used with three different regimes of spectral resolving power, R ≡ λ/Δλ. All of the observations used in the present study were obtained using the ‘medium’ slicer scale and low-resolution grating (BL), providing integral field spectra over a contiguous field of view (FoV) of 20|${_{.}^{\prime\prime}}$|3 × 16|${_{.}^{\prime\prime}}$|5 covering the common wavelength range of 3530–5530 Å with resolving power 〈R〉 = 1800.

Given the relatively small solid angle of the KCWI FoV, and the total integration time desired for each pointing of ∼5 h, it was necessary to choose the KCWI pointings carefully. In general, we chose KCWI pointings to maximize the number of previously identified KBSS galaxies with 2 ≲ z ≲ 3.4 within the FoV, so that the KCWI spectra would include the |$\rm Ly\alpha$| line as well as many other rest-frame far-UV transitions. Most of the targeted galaxies within each pointing were observed as part of KBSS in both the optical (Keck/LRIS; Oke et al. 1995; Steidel et al. 2004) and the near-IR (Keck/MOSFIRE; McLean et al. 2012; Steidel et al. 2014). The pointings were chosen so that the total sample of KBSS catalogue galaxies observed would span the full range of galaxy properties represented in the KBSS survey in terms of stellar mass (M_*), star formation rate (SFR), |$\rm Ly\alpha$| emission strength, and rest-optical nebular properties. Most of the KCWI pointings were also directed within the regions of the KBSS survey fields that have been observed at high spatial resolution by HST.

As of this writing, the KBSS–KCWI survey comprises 39 pointings of KCWI, including observations of 101 KBSS galaxies, of which 91 have 2 ≤ z ≤ 3.4 placing |$\rm Ly\alpha$| within the KCWI wavelength range. In this work, we focus only on galaxies without obvious spectroscopic or photometric evidence for the presence of an AGN, and have therefore excluded 14 objects with spectroscopic evidence for the presence of AGN, to be discussed elsewhere. The remaining objects show no sign of significant AGN activity – e.g. they lack emission lines of high ionization species in the rest-UV KCWI and existing LRIS spectra, their nebular line ratios in the rest-frame optical are consistent with stellar excitation based on spectra taken with Keck/MOSFIRE (see e.g. Steidel et al. 2014), they lack power-law SEDs in the near-IR–mid-IR, etc. Because measurements of galaxy morphology are important to the analysis, we considered only the subset of the star-forming (non-AGN) galaxies that have also been observed at high spatial resolution using HST imaging or Keck/OSIRIS IFU spectroscopy behind adaptive optics (Section 3.2).

In addition to the known KBSS targets, many of the KCWI pointings include continuum-detected serendipitous galaxies whose KCWI spectra constitute the first identification of their redshifts. A total of 50 galaxies with z > 2 (most of which are fainter in the optical continuum than the KBSS limit of R = 25.5) have been identified. We have included 10 such objects in our analysis sample, based on their having HST observations with sufficient S/N for determination of morphology.

A minimum total integration of 2.5 h (5 h is more typical) at the position of a galaxy in the KCWI data cube was also imposed, in order to ensure the relative uniformity of the data set. Ultimately, after inspection of the high-resolution images (Section 3.2), six galaxies were removed because of source ambiguity or obvious contamination from nearby unrelated objects in the images, and two were excluded because they were not sufficiently resolved by HST to measure their position angle (PA) reliably (see Section 4).

The final sample to be considered in this work contains 59 galaxies, listed in Table 1. The redshifts given in Table 1 are based on MOSFIRE nebular spectra for 38 of the 59 galaxies, which have a precision of ∼±20 km s|$^{-1}\,$| and should accurately reflect the galaxy systemic redshift. In the remaining cases, features in the rest-frame UV spectra including |$\rm Ly\alpha$| emission and strong interstellar absorption features (e.g. Si ii, Si iv, C ii, C iv, O i) were used to estimate the systemic redshift of the galaxy, using the calibration described by Chen et al. (2020b). Briefly, the method uses the statistics – based on several hundred KBSS galaxies with both nebular and UV observations – of the velocity offsets between nebular redshifts and redshifts defined by UV spectral features for samples divided by their UV spectral morphology, i.e. |$\rm Ly\alpha$| emission only, |$\rm Ly\alpha$| emission and interstellar absorption, or interstellar absorption only. The mean offsets for the appropriate sub-sample were applied to the UV-based redshifts in cases where nebular redshifts are not available; systemic redshifts obtained using such calibrations have an uncertainty of ≃ 100 km s|$^{-1}\,$| when the rest-UV spectra are of high quality (see e.g. Steidel et al. 2018). Fig. 1 shows the redshift distribution of the KCWI sample, which has z_med = 2.29 ± 0.40 (median and standard deviation), for which the conversion between angular and physical scales is 8.21 pkpc/″ with our assumed cosmology.

Reliable SED fits are available for 56 of the 59 galaxies using the BPASSv2.2 stellar population synthesis model (Stanway & Eldridge 2018) and SMC extinction curve. This choice of SED model has been shown to predict internally consistent stellar mass (M_*) and SFR for high-redshift galaxies having properties similar to those in our sample (see e.g. Steidel et al. 2016; Strom et al. 2017; Theios et al. 2019), i.e. 8.5 ≲ log (M_*/M_⊙) ≲ 11 and 1 ≲ SFR/(M_⊙ yr⁻¹) ≲ 100. The distributions of M_* and SFR (Fig. 2) are similar to those of the full KBSS galaxy sample, albeit with a slight over-representation of log (M_*/M_⊙) ≲9 galaxies¹.

In order to facilitate comparison of the |$\rm Ly\alpha$| emission line strength of sample galaxies with those in the literature, Table 1 includes the rest-frame |$\rm Ly\alpha$| equivalent width (⁠|$W_{\lambda }(\rm Ly\alpha)$|⁠) based on extraction of 1D spectra from the KCWI data cubes over a spatial aperture defined by the extent of the UV continuum light of each galaxy.² The values of |$W_{\lambda }(\rm Ly\alpha)$| in Table 1 were measured using the method described in Kornei et al. (2010) (see also Reddy & Steidel 2009), where positive values indicate net emission and negative values indicate net absorption. Aside from a slight over-representation of galaxies with the strongest |$\rm Ly\alpha$| emission (⁠|$W_{\lambda }(\rm Ly\alpha)$|≳ 40 Å), the sample in Table 1 is otherwise typical of UV-continuum-selected galaxies in KBSS by construction. The |$W_{\lambda }(\rm Ly\alpha)$| distribution for the sample in Table 1 is shown in Fig. 3. In addition to |$W_{\lambda }(\rm Ly\alpha)$|⁠, we also measured the total |$\rm Ly\alpha$| flux (F_Ly α) and the ratio between the blue- and redshifted components of emission for the entire |$\rm Ly\alpha$| halo (F_Ly α(blue)/F_Ly α(red)), which are discussed further in Section 6.

3 OBSERVATIONS AND REDUCTIONS

3.1 KCWI

The KCWI data discussed in the present work were obtained between 2017 September and 2020 November, in all cases using the medium-scale slicer made up of 24 slices of width 0|${_{.}^{\prime\prime}}$|69 and length 20|${_{.}^{\prime\prime}}$|3 on the sky. The instrumental setup uses the BL volume phase holographic grating with an angle of incidence that optimizes the diffraction efficiency near 4200 Å, with the camera articulation angle set to record the spectra of each slice with a central wavelength of ∼4500 Å. A band-limiting filter was used to suppress wavelengths outside of the range of 3500–5600 Å to reduce scattered light; the useful common wavelength range recorded for all 24 slices in this mode is 3530–5530 Å, with a spectral resolving power ranging from R ≃ 1400 at 3530 Å to R ≃ 2200 at 5530 Å. At the mean redshift of the sample (〈z〉 = 2.42), |$\rm Ly\alpha$| falls at an observed wavelength of ∼4160 Å, where R ∼ 1650.

The E2V 4k × 4k detector was binned 2 × 2, which provides spatial sampling along slices of 0|${_{.}^{\prime\prime}}$|29 pix⁻¹. Because each slice samples 0|${_{.}^{\prime\prime}}$|68 in the dispersion direction, the effective spatial-resolution element is rectangular on the sky, with an aspect ratio of ∼2.3: 1. We adopted the following approach to the observations, designed to ensure that the slicer geometry with respect to the sky is unique on each 1200-s exposure so that the effective spatial resolution on the final stacked data cube is close to isotropic. Typically, a total integration of ∼5 h is obtained as a sequence of 15 exposures of 1200 s, each obtained with the sky PA of the instrument rotated by 10–90 degrees with respect to adjacent exposures. Each rotation of the instrument FoV is accompanied by a small offset of the telescope pointing before the guide star is reacquired. In this way, a given sky position is sampled in 15 different ways by the slicer.

3.1.1 KCWI data reduction

Each 1200-s exposure with KCWI was initially reduced using the data reduction pipeline (DRP) maintained by the instrument team and available via the Keck Observatory web site.³ The DRP assembles the 2D spectra of all slices into a 3D data cube (with spaxels of 0|${_{.}^{\prime\prime}}$|29 × 0|${_{.}^{\prime\prime}}$|69, the native scale) using a suite of procedures that can be customized to suit particular applications. The procedures include cosmic ray removal, overscan subtraction and scattered light subtraction, wavelength calibration, flat-fielding, sky subtraction, differential atmospheric refraction (DAR) correction, and flux calibration. Wavelength calibration was achieved using ThAr arc spectra obtained using the internal calibration system during the afternoon prior to each observing night. Flat-fielding was accomplished using spectra of the twilight sky at the beginning or end of each night, after dividing by a b-spline model of the solar spectrum calculated using the information from all slices. For each frame, the sky background was subtracted using the sky-modelling feature in the DRP, after which the sky-subtracted image (still in the 2D format) is examined in order to mask pixels, in all 24 slices, that contain significant light from sources in the field. The frame was then used to make a new 2D sky model using only unmasked pixels, and the sky-subtracted image is reassembled into a wavelength-calibrated (rebinned to 1 Å per wavelength bin) data cube, at which time a variance cube, an exposure cube, and a mask cube are also produced.

Next, we removed any remaining low-frequency residuals from imperfect sky background subtraction by forming a median-filtered cube after masking obvious continuum and extended emission line sources using a running 3D boxcar filter. The typical dimensions of the filter are 100 Å (100 pixels) in the wavelength direction, 16 pixels (4|${_{.}^{\prime\prime}}$|6) along slices, and 1 pixel (0|${_{.}^{\prime\prime}}$|69) perpendicular to the slices, with the last ensuring slice-to-slice independence. Minor adjustments to the filter dimensions were made as needed. If the running boxcar encounters a large region with too few unmasked pixels to allow a reliable median determination, then the pixel values in the filtered cube were interpolated from the nearest adjacent regions for which the median was well determined. Finally, the median-filtered cube for each observed frame was subtracted from the data. We found that this method proved effective for removing scattered light along slices caused by bright objects within the KCWI FoV.

Because neither Keck II nor KCWI has an atmospheric dispersion corrector, each cube is corrected for DAR (i.e. apparent position of an object as a function of wavelength) using the elevation and parallactic angle at the midpoint of the exposure and a model of the atmosphere above Maunakea. Finally, each cube was flux-calibrated using observations obtained with the same instrument configuration of one or more spectrophotometric standard stars selected from a list recommended by the KCWI documentation.

Prior to stacking reduced data cubes of individual exposures covering the same sky region, they must be aligned spatially and rotated to account for differences in the PA of the instrument with respect to the sky on different exposures. To accomplish this, we averaged the DAR-corrected cubes along the wavelength axis to create pseudo-white-light images, rotated each to the nominal sky orientation (N up and E left) based on the World Coordinate System (WCS) recorded in the header, and cross-correlated the results to determine the relative offsets in RA and Dec, which we found to have a precision of ≃0.03 arcsec (root-mean-square, RMS). The offsets were applied to the WCS in the header of each cube, and all cubes for a given pointing were then resampled to a common spatial grid with spatial sampling of 0|${_{.}^{\prime\prime}}$|3 × 0|${_{.}^{\prime\prime}}$|3 using a 2D drizzle algorithm in the Montage package,⁴ with drizzle factor of 0.7. If the wavelength grid is different among cubes, the spectrum in each spaxel was resampled to a common grid using cubic spline. Finally, we stacked the resampled cubes by averaging, weighted by exposure time. A white light image of the final stacked data cube was used to determine small corrections to the fiducial RA and Dec to align with other multiwavelength data of the same region.

3.2 High spatial-resolution imaging

For galaxies in the mass range of our sample at z ∼ 2.5, the typical half-light diameter is ∼4 pkpc, which corresponds to ≃0|${_{.}^{\prime\prime}}$|5 (Law et al. 2012b). To obtain reliable measurements of the orientation of the galaxy’s projected major and minor axes (see Section 4), high-resolution images are crucial. We gathered existing data from two sources: space-based optical or near-IR images from the HST/ACS or the HST/WFC3; alternatively, H α maps obtained using the Keck/OSIRIS integral field spectrometer (Larkin et al. 2006) behind the WMKO laser guide star adaptive optics system, which typically provides spatial resolution of 0|${_{.}^{\prime\prime}}$|11–0|${_{.}^{\prime\prime}}$|15 (Law et al. 2007, 2009, 2012a).

The HST/WFC3-IR images in this work were obtained using either the F140W or F160W filters, from programs listed in Table 2, with spatial resolution of ∼0|${_{.}^{\prime\prime}}$|16 and∼0|${_{.}^{\prime\prime}}$|18, respectively. The HST/ACS images were taken in either the F606W or F814W filters, with spatial resolution of ∼0|${_{.}^{\prime\prime}}$|09 (full width at half-maximum; FWHM). In all cases, overlapping exposures were aligned and combined using DrizzlePac.⁵ For the Q2343 field, where galaxies have been observed with comparably deep observations in two filters, we selected the image with the smaller estimated uncertainty in the measured PA (see Section 4.2). The Keck/OSIRIS |$\rm H\alpha$| maps of six galaxies included in the sample are presented by Law et al. (2007), Law et al. (2009), Law et al. (2012a), and Law et al. (2018); details of the observations and data reduction can be found in those references.

4 GALAXY AZIMUTHAL ANGLE

4.1 Motivation

The physical state of the CGM is controlled by the competition between accretion of new material on to the galaxy and outflows driven by processes originating at small galactocentric radii. In what one might call the ‘classic’ picture (see Veilleux, Cecil & Bland-Hawthorn 2005; Tumlinson et al. 2017 for review articles), based on observations of nearby starburst galaxies, inflowing and outflowing gas occurs preferentially along the major and minor axes, respectively, as in the schematic diagram in Fig. 4. It is well established around star-forming galaxies at z ≲ 1 that outflows driven by energy or momentum from stellar feedback (radiation pressure from massive stars, supernovae) originating near the galaxy centre tend to escape along the direction that minimizes the thickness of the ambient interstellar medium (ISM) that would tend to slow or prevent outflows from escaping the galactic disc. Meanwhile, accretion of cool gas from the intergalactic medium (IGM) is believed to occur in quasi-collimated filamentary flows that carry significant angular momentum and thus would tend to approach the inner galaxy in a direction parallel to the disc plane (Nelson et al. 2019; Péroux et al. 2020). If feedback processes are sufficiently vigorous, the accretion is also most likely along directions that do not encounter strong outflows. When projected on to the plane of the sky, these considerations would tend to intersect outflowing material along the polar direction (the projected minor axis) and accreting material when the azimuthal angle (see Fig. 4) is closer to the PA of the projected major axis.

At redshifts z ≲ 1, there is strong support for this general geometric picture from the statistics of the incidence and strength of rest-UV absorption lines observed in the spectra of background sources, as a function of the azimuthal angle of the vector connecting the sightline and the centre of the galaxy. As might be expected from a picture similar to Fig. 4, sightlines with ϕ close to the minor axis intersect gas with a large range of velocities, which kinematically broadens the observed absorption complexes comprising many saturated components (e.g. Bordoloi et al. 2011; Bouché et al. 2012; Kacprzak et al. 2012; Schroetter et al. 2019). For azimuthal angles ϕ close to the major axis, the absorption features are strong due to the high covering fraction and column density of low ionization gas near the disc plane and broadened by the kinematics of differential rotation. The same picture has been supported by cosmological simulations where the gas metallicity, radial velocity, and the amount of outflowing mass all show significant differences along the galaxy major and minor axes (Nelson et al. 2019; Péroux et al. 2020).

It remains unclear, however, whether this geometric picture should be applied to galaxies at z ∼ 2, where a large fraction of galaxies, particularly those with log(M_*/M_⊙) ≲ 10, appear to be ‘dispersion dominated’, where the rotational component of dynamical support (V_rot) is significantly smaller than that of apparently random motion (σ) (Förster Schreiber et al. 2009; Law et al. 2009), and where the central dynamical mass of the galaxy may be dominated by cold gas rather than stars – and therefore any disc would be highly unstable. For such galaxies, there is not always a clear connection between the morphology of starlight and the principal kinematic axes (e.g. Erb et al. 2004; Law et al. 2012b). Meanwhile, a WFC3 survey conducted by Law et al. (2012b) for similar galaxies strongly supports 3D triaxial morphology, instead of inclined disc, to be the most suitable model to describe these galaxies. The fact that the vast majority of galaxies with DTB spectra at z > 2 have systematically blueshifted interstellar absorption lines and systematically redshifted |$\rm Ly\alpha$| emission suggests that outflows cannot be confined to a small range of azimuth (Shapley et al. 2003; Steidel et al. 2010; Jones, Stark & Ellis 2012 etc.).

The spatial and spectral distribution of the |$\rm Ly\alpha$| emission surrounding galaxies may provide crucial insight into the degree of axisymmetry and the dominant magnitude and direction of gas flows in the CGM. In any case, IFU observations of |$\rm Ly\alpha$|⁠, where both the geometry and the kinematics of the extended emission can be mapped, complement information available from absorption line studies.

4.2 Methods

In order to obtain reliable measurements of ϕ, it is important to measure |$\rm PA_0$| accurately and consistently. We used up to four different methods to determine |$\rm PA_0$| for each galaxy. The choice of method depends on the information available; these are briefly summarized as follows:

• Sérsic profile fitting of HST images: In this method, we fit a 2D Sérsic profile (Sérsic 1963) to the host galaxy by using galfit (Peng et al. 2002, 2010) and determined |$\rm PA_0$| from the best-fitting model parameters. The point spread function (PSF) was measured by selecting stellar sources over the full HST pointing using the star classifier in SExtractor, which calculates a ‘stellarity index’ for each object based on a neural network. We then examined sources with the highest 3 per cent stellarity indices by eye and normalized and stacked them to form an empirical PSF. Our fiducial model consists of a 2D elliptical Sérsic profile convolved with the PSF. We also included the first-order Fourier mode to handle the asymmetric morphology in most cases. However, over- or under-fitting can cause failure of convergence or unreasonably large fitting errors and residuals. In most cases, the cause of the failure and the required adjustment are obvious in the original galaxy image and the model residual. For example, if the residual reveals an additional source, we would add an additional Sérsic component or a simple scaled PSF to the model, depending on the size of the additional source. Meanwhile, if the primary source shows a triangular morphology, we would add the third-order Fourier mode⁶ associated with the Sérsic profile. However, in certain cases, obtaining a successful fit requires experimenting with the model by adding or removing certain degrees of freedom. The rule of thumb is that we add or remove degrees of freedom one at a time and adopt the adjustment if it makes the fit converge or significantly diminishes the reduced χ² and the fitted error of |$\rm PA_0$|⁠. In the end, 23 galaxies were fit with the fiducial model. Nine galaxies (Q0100-C7, Q0821-MD36, Q1009-BX222, Q1009-D15, Q1623-BX453-CS3, Q1700-BX729-CS9, Q2343-BX418, Q2343-BX418-CS8, and Q2343-BX660) were fit without the Fourier modes. Five galaxies (Q1009-BX218, Q1700-BX490, Q1700-BX710, Q1700-BX729, and Q2343-BX513-CS7) were fit with additional sources. Three galaxies (Q2343-BX391, Q2343-BX587-CS3, and Q2343-BX587-CS4) were fit with third-order Fourier modes. Galaxies with unsuccessful fits require using alternative methods, detailed below. Successful Sérsic fits were obtained for 40 of 59 galaxies, shown in Figs 5 and 7; galaxies with successful applications of this method tend to be isolated and to have a dominant central high SB component.

• Second moment of pixel intensity on HST images: This method calculates the flux-weighted pixel PA arithmetically derived from the second moment of the pixel intensity,

  $$\begin{eqnarray*} \tan (2 {\rm PA_0}) = \frac{2 \langle xy \rangle }{\langle x^2 \rangle - \langle y^2 \rangle }, \end{eqnarray*}$$

  (2)

  where x and y are the pixel positions relative to the centre of the galaxy in the x and y directions, and ‘〈...〉’ indicates the arithmetic mean of the pixel coordinates weighted by pixel flux. Galaxy centres were defined as the point where 〈x〉 = 0 and 〈y〉 = 0. Especially for the galaxies that are morphologically complex, we found that the |$\rm PA_0$| measurements using this method are sensitive to the SB threshold used to define the outer isophotes and to the spatial resolution of the HST image. Therefore, all ACS images (which have higher spatial resolution than those taken with WFC3-IR) were convolved with a 2D Gaussian kernel to match the PSF to those of the WFC3-IR/F160W images. Among 59 galaxies, |$\mathrm{PA}_{0}$| of 58 galaxies can be otained using this method. They are shown in Figs 5, 6, and 7. We tested various SB thresholds and found that a threshold of 50 per cent of the peak SB after convolution provides the most consistent measurements between the |$\rm PA_0$| measured using (i) and (ii) has an rms of 10.8 degrees.

• Second moment of pixel intensity of OSIRIS H α maps: This method uses the same algorithm as in (ii), applied to Keck/OSIRIS H α maps rather than HST continuum images. The SB threshold follows Law et al. (2009). There are six galaxies whose major axes were determined using this method, including one (Q1623-BX453) whose F160W image is unresolved. For the remaining five galaxies, the rms between |$\rm PA_0$| measured from the HST images and this method is ∼15 degrees. All six galaxies are shown in Fig. 7.

• H α Kinematics: As the only galaxy in this sample that demonstrates not only rotational kinematics from H α emission but also clear disc morphology (Law et al. 2012a), the galaxy major axis of Q2343-BX442 is defined to be perpendicular to its rotational axis (Fig. 8). Because of its complex morphology, we did not apply method (i) to this galaxy. Method (ii) is likely dominated by the inner spiral structure and leads to a |$\rm PA_0$| that is ∼90 degrees apart from that determined from kinematics, while method (iii) is consistent with this method within 15 degrees.

Table 1 lists the adopted measurement of |$\rm PA_0$| for each galaxy in the sample, as well as the methods used to determine it. In cases where multiple methods were applied, we used the average |$\rm PA_0$| of the OSIRIS and HST measurements, where the latter value is an average of the results obtained using methods (i) and (ii). This way, more weight is given to method (iii) results when available, since it utilizes entirely independent data from a different instrument. The robustness of |$\rm PA_0$| measurements is discussed in the following two subsections.

4.3 Relationship between the kinematic and morphological major axes

An important issue for the interpretation of morphological measurement of |$\rm PA_0$| is the extent to which it is likely to be a proxy for the kinematic major axis. Two of the 17 galaxies in Fig. 6 (Q2343-BX389 and Q2343-BX610) were studied as part of the SINS IFU survey by Förster Schreiber et al. (2018, hereafter FS18) and were subsequently re-observed with the benefit of adaptive optics by FS18; the latter found that both have clear rotational signatures and that the inferred differences between their morphological major axis |$\rm PA_0$| and the kinematic major axis |$\rm PA_{\rm kin}$| are |$\Delta \rm PA\equiv |\rm PA_0 - \rm PA_{\rm kin}| = 2\pm 5$| deg and 27 ± 10 deg for Q2343-BX389 and Q2343-BX610, respectively. From their full sample of 38 galaxies observed at AO resolution, FS18 found 〈ΔPA〉 = 23 degrees (mean), with median ΔPA_med = 13 degrees.

Among the six galaxies that also have H α velocity maps from Law et al. (2009) (used for method (iii) above), two (Q0449-BX93 and Q1623-BX453) show no clear rotational signature, so that |$\rm PA_{{\rm kin}}$| is indeterminate. Small amounts of velocity shear were observed for Q2343-BX513 and Q2343-BX660, both having |$\rm PA_{\rm kin}$| consistent with |$\rm PA_0$|⁠.⁷ However, for Q2343-BX418, the implied |$\rm PA_{\rm kin}$| is nearly perpendicular to its morphology-based |$\rm PA_0$|⁠. For Q1700-BX490, the kinematic structure is complicated by the presence of two distinct components: the brighter, western component appears to have |$\rm PA_{\rm kin} \simeq 20$| degrees, which would be consistent with |$\rm PA_0$| measured from its |$\rm H\alpha$| intensity map; however, if the eastern component, which has a slightly blueshifted velocity of |$\sim 100 \rm km~s^{-1}\,$|⁠, is included as part of the same galaxy, we found |$\rm PA_0 = 86$| degrees.

As summarized in Table 1 and Figs 5–7, most of the galaxy sample has |$\rm PA_0$| measured using two or more of the methods described above, and generally the morphologically determined major axis |$\rm PA_0$| agrees with one another to within ∼10 degrees. We caution that no high spatial-resolution kinematic information is available for most of the sample; based on the subset of nine that do have such measurements, approximately two-thirds show reasonable agreement between |$\rm PA_0$| and |$\rm PA_{\rm kin}$|⁠.

4.4 Distribution and robustness of |$\rm PA_0$| measurements

To estimate the systematic uncertainty of the |$\rm PA_0$| measurements, we compare values measured using different methods for the same objects in Fig. 9. Between methods (i) and (ii), which are both based on HST images, the measured |$\rm PA_0$| for 40 galaxies are well centred along the 1:1 ratio, with rms ≃ 10.8 degrees. Values of |$\rm PA0$| measured using method (iii) are based on spectral line maps from an IFU rather than continuum light in a direct image, so that the 5 |$\rm PA_0$| values measured using methods (i) and (iii) exhibit larger scatter, with rms ≃15.6 degrees relative to a 1:1 ratio. Therefore, we conclude that the uncertainty in final |$\rm PA_0$| measurements is ≲15 degrees.

Fig. 10 shows the normalized kernel density estimator (KDE) and the individual measurements of |$\rm PA_0$|⁠. There is an apparent excess in the occurrence rate of values between |$\rm PA_0 \sim 10-40^\circ$|⁠. To evaluate its possible significance, we conducted 1000 Monte Carlo realizations of a sample of 58 galaxies with randomly assigned |$\rm PA_0$|⁠; we find that there is an ≃5 per cent probability that a similar excess is caused by chance, thus is not statistically significant, and is consistent with that expected from sample variance.

Furthermore, because the HST and OSIRIS images used to measure |$\rm PA_0$| were rotated to the nominal North up and East left orientation prior to measurement, we tested whether significant bias might result from the choice of pixel grid by re-sampling all images with pixel grids oriented in five different directions. We found that the values measured for |$\rm PA_0$| were consistent to within 10 degrees (rms).

In order to characterize the PSF of our final KCWI data cubes, and to show that it is effectively axisymmetric, Fig. 11(b) shows various composite images of the galaxy sample. Prior to forming the stack, each individual data cube was rotated to align the galaxy major (minor) axis with the X (Y) coordinate of the stacked image. The KCWI stacked galaxy continuum image [panel (i)] was made by integrating each aligned data cube along the wavelength axis over the range 1224 ≤ λ₀/Å ≤ 1236 in the galaxy rest frame (i.e. |$2000 \lt \Delta v/\rm km~s^{-1}\, \le 5000$|⁠). This integration window was chosen to be representative of the UV continuum near |$\rm Ly\alpha$| without including the |$\rm Ly\alpha$| line itself and to be unaffected by |$\rm Ly\alpha$| absorption from the IGM. The centre of the galaxy is determined with high confidence (see Section 5.1) by fitting a 2D Gaussian function to the individual galaxy continuum image. A similar approach – aligning the principal axes in the high-resolution HST images prior to stacking – was used to produce the HST stacked continuum image shown in panel (ii), except that the centres in the HST images were measured using the first moment of pixel intensity. Both stacks were conducted in units of observed SB. The FWHM of the stacked HST image along the major and minor axes are 0|${_{.}^{\prime\prime}}$|55 and 0|${_{.}^{\prime\prime}}$|35, respectively. The stacked KCWI continuum image is well reproduced by the convolution of the |$\rm PA_0$|-aligned stacked HST image [panel (ii) of Fig. 11(b)] with an axisymmetric 2D Gaussian profile with FWHM = 1|${_{.}^{\prime\prime}}$|02 [see panels (iii) and (iv) of Fig. 11(b)].

Comparing the convolution of HST and KCWI [panel (iii)] with the best direct fit of a symmetric Gaussian profile to the KCWI continuum image (⁠|$\rm FWHM = 1{{_{.}^{\prime\prime}} }21$|⁠) [panel (v)] shows that even at the ≃1|${_{.}^{\prime\prime}}$|02 resolution of KCWI one can clearly distinguish the major axis elongation. The residual map assuming a symmetric Gaussian profile [panel (vi)] shows a clear dipole residual compared to panel (iv). Thus, Fig. 11 shows that (1) the PSF of the KCWI cubes is axisymmetric and thus has not introduced a bias to the azimuthal light distribution measurements and (2) the spatial resolution is sufficient to recognize non-axisymmetry even on sub-arcsec angular scales of the continuum light.

5 ANALYSES

5.1 |$\rm Ly\alpha$| spatial profile

To study the dependence of the |$\rm Ly\alpha$| emission profile on galaxy azimuthal angle, we first analyze the |$\rm Ly\alpha$| SB profile as a function of the impact parameter (or transverse distance, D_tran). Fig. 12 compares the stacked continuum and continuum-subtracted narrow-band |$\rm Ly\alpha$| emission, composed in the same way as in Fig. 11. The integration window of the |$\rm Ly\alpha$| image is −700 < Δv/(kms⁻¹) ≤ 1000 (1213–1220 Å) to include most of the |$\rm Ly\alpha$| emission (as shown later in Section 5.2). Continuum subtraction throughout this work was done spaxel-by-spaxel in the data cube of each target galaxy by subtracting the average flux density in two windows flanking the position of |$\rm Ly\alpha$|⁠, with 2000 < |cΔv_sys/km s⁻¹| < 5000, where Δv_sys is the velocity separation relative to rest-frame |$\rm Ly\alpha$| (i.e. two windows each of width ≃ 12 Å in the rest frame [1195–1207] and [1224–1236] Å. Similar to Section 4.4, the centre of the galaxy in the KCWI data was determined by fitting a 2D Gaussian profile to the KCWI continuum image. The fitting box was chosen by eye to include all of the signal from the galaxy but excluding nearby contamination. Despite the arbitrary box size, we found that the derived centroid is very robust: varying the box size by 4 pixels (1|${_{.}^{\prime\prime}}$|2), the fit result does not change by more than 0|${_{.}^{\prime\prime}}$|01. The typical fitting error propagated from the reduced χ² is also ∼0|${_{.}^{\prime\prime}}$|01 (median), i.e. much smaller than the seeing disc.

A 2D Gaussian fit to the profiles finds that the FWHM values are 1|${_{.}^{\prime\prime}}$|279 ± 0|${_{.}^{\prime\prime}}$|003 (major axis) × 1|${_{.}^{\prime\prime}}$|182 ± 0|${_{.}^{\prime\prime}}$|003 (minor axis) (or an |$\sim 8{{\ \rm per\ cent}}$| difference) for the continuum emission and 1|${_{.}^{\prime\prime}}$|689 ± 0|${_{.}^{\prime\prime}}$|005 × 1|${_{.}^{\prime\prime}}$|705 ± 0|${_{.}^{\prime\prime}}$|005 (⁠|$\lt 1{{\ \rm per\ cent}}$| difference between major and minor axes) for the |$\rm Ly\alpha$| emission. Therefore, the |$\rm Ly\alpha$| emission in the stacked image is both more symmetric and more spatially extended than the continuum emission.

Fig. 13 shows the median |$\rm Ly\alpha$| SB as a function of D_tran (red). Each point represents the median SB of a bin of pixels with ΔD_tran = 0.1 pkpc. The |$\rm Ly\alpha$| SB profile falls off much more slowly than the continuum (grey). Following Wisotzki et al. (2016), we fit the |$\rm Ly\alpha$| SB profile with a two-component model – a compact ‘core’ component and an extended ‘halo’ component. Both components are exponential profiles convolved with the KCWI PSF with the amplitudes, exponential radii, and a uniform background term as free parameters. Further details of the model fitting will be described in a future work. The best-fitting exponential radii |$r_\mathrm{exp} = 3.71^{+0.06}_{-0.04}$| pkpc and |$15.6^{+0.5}_{-0.4}$| pkpc. The r_exp of the halo component is close to Steidel et al. (2011) (for KBSS galaxies observed with narrow-band |$\rm Ly\alpha$| imaging), which found that the median-stacked |$\rm Ly\alpha$| profile has r_exp = 17.5 pkpc but slightly more extended than Wisotzki et al. (2018) for SF galaxies at z > 3 (see also Matsuda et al. 2012; Momose et al. 2014; Leclercq et al. 2017). Dividing the SB profiles into two subsamples with 0° ≤ ϕ < 45° (purple) and 45° ≤ ϕ ≤ 90° (orange) that represent the galaxy major and major axes, respectively, one can see that the resulting profiles are consistent with one another to within 1σ, or within |$\lesssim 2\times 10^{-19}\, \mathrm{~erg~s}^{-1}\, \mathrm{cm}^{-2}\, \mathrm{arcsec}^{-2}$|⁠. The possible exception is at the smallest projected distances (D_tran < 1 arcsec, or ≲8 pkpc), where the |$\rm Ly\alpha$| emission is marginally enhanced; if real, the difference in profiles (the asymmetry) represents |$\lt 2{{\ \rm per\ cent}}$| of the total |$\rm Ly\alpha$| flux. Thus, the composite |$\rm Ly\alpha$| intensity is remarkably symmetric, suggesting an overall lack of a strong statistical connection between the morphology of the starlight and that of the extended |$\rm Ly\alpha$| emission surrounding individual star-forming galaxies at z ∼ 2–3.

5.2 Cylindrical projection of 2D spectra

The similarity of the |$\rm Ly\alpha$| SB profile along galaxy major and minor axes suggests that extended |$\rm Ly\alpha$| emission depends little on the galaxy orientation. However, the KCWI data cubes allow for potentially finer discrimination through examination of both the SB and the kinematics of |$\rm Ly\alpha$| emission as a function of projected galactocentric distance. To facilitate such comparison, we introduce ‘cylindrical projections’ of 2D |$\rm Ly\alpha$| emission. The basic idea behind cylindrical projection, illustrated in Fig. 14, is to provide an intuitive visualization of spatial and spectral information simultaneously.

Compared to the standard 2D spectrum one obtains from slit spectroscopy, the cylindrical 2D spectrum replaces the 1D spatial axis (i.e. distance along a slit) with projected distance by averaging spaxels in bins of D_tran or, equivalently, θ_tran. When projected as in the right-hand panel of Fig. 14, it can also be viewed as the |$\rm Ly\alpha$| spectrum at every projected radial distance (averaged, in this case, over all azimuthal angles) or as the average radial profile at each slice of wavelength or velocity.

Fig. 14 shows the stacked cylindrical 2D spectrum formed by averaging the continuum subtracted data cubes at wavelengths near rest-frame |$\rm Ly\alpha$| for all 59 galaxies in Table 1. This composite cylindrical 2D spectrum, analogous to a ‘down-the-barrel’ |$\rm Ly\alpha$| spectrum in 1D, but evaluated as a function of galactocentric distance, shows that the |$\rm Ly\alpha$| emission line comprises distinct redshifted and blueshifted components extending to ±1000 km s|$^{-1}\,$| with respect to v_sys = 0, with a minimum close to v_sys = 0. The vast majority of individual galaxies, and therefore also the average in the stacked profile, has F_Ly α(blue)/F_Ly α(red) ≃ 0.3 and is thus ‘red peak dominated’. The two-component spectral morphology extends to at least θ_tran ≃ 3 arcsec or D_tran ≃ 25 pkpc. This overall spectral morphology is most readily explained as |$\rm Ly\alpha$| photons being resonantly scattered by outflowing material, whereby redshifted photons scattered from the receding (opposite) side are more likely to escape in the observer’s direction than blueshifted photons (e.g. Pettini et al. 2001; Steidel et al. 2010; Dijkstra 2014).

As D_tran increases and the |$\rm Ly\alpha$| SB decreases exponentially, the two |$\rm Ly\alpha$| peaks become less distinct and merge into a symmetric ‘halo’ centred on Δv = 0. The vast majority of the |$\rm Ly\alpha$| emission is within −700 < Δv/km s⁻¹ < 1000. However, we caution that the apparent blue edge at Δv ∼ −700 km s⁻¹ of |$\rm Ly\alpha$| emission in this composite 2D spectrum is likely caused by continuum over-subtraction resulting from the relatively simple technique that we used in this work. The continuum subtraction assumed a linear interpolation of the continuum spectrum underneath the |$\rm Ly\alpha$| emission (see Section 5.1), which tends to over-estimate the continuum flux blueward of the systemic redshift of the galaxy due to intrinsic |$\rm Ly\alpha$| absorption in the stellar spectrum and residual effects of the often-strong |$\rm Ly\alpha$| absorption damping wings on which |$\rm Ly\alpha$| emission is superposed. To improve on this would require a more sophisticated continuum subtraction method in the inner ≃1|${_{.}^{\prime\prime}}$|0 of the galaxy profile; however, since most of the remainder of this work will involve comparison of 2D cylindrical projections with one another, the imperfections in continuum subtraction at small θ_tran are unlikely to affect the results.

5.3 Dependence on azimuthal angle of cylindrically projected 2D (CP2D) spectra

To investigate how |$\rm Ly\alpha$| emission depends on the galaxy azimuthal angle, we split each CP2D spectrum of individual galaxies averaged over two independent bins of azimuthal angle (ϕ) with respect to the galaxy’s major axis: 0° ≤ ϕ < 45° (‘major axis’) and 45° ≤ ϕ ≤ 90° (‘minor axis’), as in Section 5.1 and Fig. 13. The CP2D stacks covering these azimuth ranges were combined separately to form CP2D composites that we refer to as ‘major axis’ and ‘minor axis’. To reveal subtle differences in SB and/or velocity along these two directions, we subtracted one from the other – Fig. 15 shows the result.

The difference between the CP2D spectra along the major and minor axes (top right of Fig. 15) shows excess emission along the galaxy major axis at Δv ≃ +300 kms⁻¹ – consistent with the velocity of the peak of the redshifted component in the full composite CP2D spectrum – within θ_tran ≲ 2″  (D_tran ≲ 15 pkpc). The |$\rm Ly\alpha$| flux of this asymmetric component of |$\rm Ly\alpha$| amounts to ≲2 per cent of the total and has a peak intensity ≃5 per cent that of the peak of the redshifted |$\rm Ly\alpha$| component shown in Fig. 13.

The significance map in the bottom right-hand panel of Fig. 15 is based on the standard deviation of 100 independent mock CP2D stacks, each made by assigning random |$\rm PA_0$| to the galaxies in our sample before combining. The mock stacks were then used to produce a 2D map of the rms residuals evaluated in the same way as the observed data. Considering the effective resolution, the overall significance (compared to the standard deviation) of the most prominent feature in the top right-hand panel of Fig. 15 is ≃2–2.5σ per resolution element. Thus, while the residual (excess) feature may be marginally significant statistically, the level of asymmetry relative to the total |$\rm Ly\alpha$| flux is in fact very small.

5.4 The robustness of residual |$\rm Ly\alpha$| asymmetry

Despite the marginally significant detection of the excess |$\rm Ly\alpha$| emission along galaxy major axes, its robustness is subject to scrutiny. In particular, we would like to determine whether the apparent detection is typical of the population or is caused by a few outlier objects having very asymmetric |$\rm Ly\alpha$| as a function of azimuthal angle.

Fig. 16 shows histograms of the difference in integrated |$\rm Ly\alpha$| flux between the major and minor axis bins of azimuthal angle. We calculated the differences integrated over two different ranges of Δv and θ_tran; (1) the range where the excess shown in Fig. 15 is most prominent, θ_tran ≤ 2″  and |$0 \le (\Delta v)/\rm km~s^{-1}\, \le 500$|⁠, shown in the top panel, and (2) the range that encapsulates most of the redshifted component of |$\rm Ly\alpha$| θ_tran ≤ 3″  and |$0 \le (\Delta v)/\rm km~s^{-1}\, \le 1000$|⁠, shown in the bottom panel. Two galaxies – Q0142-BX165 and Q2343-BX418 – are clearly outliers in (1), while only Q0142-BX165 stands out in (2). The distribution of ΔF_Ly α for the other 56 galaxies in the sample is relatively symmetric around ΔF_Ly α = 0.

Figs 17 and 18 show the stacked |$\rm Ly\alpha$| profiles as in Fig. 15 but with the strongest outliers removed from the stack. After removing both outliers (Fig. 18), the excess |$\rm Ly\alpha$| emission along the galaxy major axis at Δv ≃ 300 km s⁻¹ becomes consistent with noise. Although Q2343-BX418 is not an extreme outlier in terms of the overall |$\rm Ly\alpha$| asymmetry of the integrated redshifted component of |$\rm Ly\alpha$| emission (bottom panel of Fig. 16), when only Q0142-BX165 is removed from the stack (Fig. 17), the composite 2D spectra still show obvious excess |$\rm Ly\alpha$| emission along the galaxy major axis, albeit with slightly reduced significance. Meanwhile, when both outliers are removed from the stack (Fig. 18), one can see an emerging excess |$\rm Ly\alpha$| emission along the minor axis, for the blue peak, with −700 ≲ Δv/(km s⁻¹) ≲ −200 and θ_tran ≲ 2|${_{.}^{\prime\prime}}$|5 with integrated significance ∼2σ. The flux of the excess blueshifted emission comprises ∼10 per cent of the total blueshifted |$\rm Ly\alpha$| flux, with a peak amplitude ∼1 per cent of the peak |$\rm Ly\alpha$| intensity (i.e. the redshifted peak).

We conducted an analysis on the blueshifted emission similar to that done for the redshifted asymmetry, with results summarized in Fig. 19. There is no obvious outlier in the difference in integrated |$\rm Ly\alpha$| flux between the major and minor axis azimuth bins except for Q0142-BX165, for which the excess again favours the major axis (i.e. it is in the direction opposite to the apparent blueshifted asymmetry identified in Fig. 18). We also consecutively removed from the |$\rm Ly\alpha$| stack galaxies with extreme excess emission along the minor axis and found no sudden and significant changes in the composite spectra. Evidently, the blueshifted excess along the minor axis, while of about the same significance as the redshifted excess in the major axis direction, is a general property of the full sample rather than a result of a small number of outliers.

In summary, we found excess emission along the galaxy major axis for the redshifted component of |$\rm Ly\alpha$| near Δv ≃ 300 km s⁻¹. However, this particular excess emission appears to be caused by galaxy outliers with extreme emission along the major axis. After removing them from the composite CP2D spectra, we found excess emission along the galaxy minor axis for the blue peak within −700 ≲ Δv/(km s⁻¹) ≲−200 that is not apparently affected by the extreme scenarios. Both detections are not particularly significant at the ∼2σ level.

5.5 A closer look at the extreme cases

Fig. 20 shows the residual between the cylindrically projected 2D spectra extracted along the major and minor axes of Q0142-BX165, as well as continuum images from both HST and KCWI. Q0142-BX165 has two comparably bright components separated by ≃ 3.5 pkpc (≃ 0|${_{.}^{\prime\prime}}$|4). Careful inspection of Keck/LRIS and Keck/MOSFIRE spectra of this system revealed no sign of an object at a different redshift. There is no significant offset between the KCWI continuum and the |$\rm Ly\alpha$| centroids (≤0.5 pix ∼0|${_{.}^{\prime\prime}}$|15 separation), indicating that the apparent directional asymmetry in the CP2D spectra is not caused by a spatial shift between the continuum and the |$\rm Ly\alpha$| emission. Instead, the narrow-band |$\rm Ly\alpha$| map shows that the |$\rm Ly\alpha$| emission is elongated approximately along the N–S direction. However, after aligning the KCWI and HST astrometry with reference to a nearby compact galaxy, we found that both the KCWI stellar continuum near |$\rm Ly\alpha$| and the narrow-band |$\rm Ly\alpha$| emission are centred near the SW component (see Fig. 20). It is possible that the SW component alone is responsible for the |$\rm Ly\alpha$| emission, in which case its |$\rm PA_0$| would be −59°,  ∼70° off from what was determined in Section 4. However, adopting |$\rm PA_0 = -59^{\circ }$| would cause BX165 to become an outlier with excess |$\rm Ly\alpha$| emission along the minor axis. Meanwhile, the elongation of the KCWI continuum aligns with the direction of the separation of the two components and is roughly consistent with the direction of the |$\rm Ly\alpha$| elongation as well. This seems to suggest that the |$\rm Ly\alpha$| elongation simply reflects the asymmetry of the continuum source, albeit on a larger angular scale; however, as shown in Fig. 6, many galaxies in the sample possess similar morphologies, but Q0142-BX165 is the only one that shows extraordinary asymmetry in |$\rm Ly\alpha$| emission. In any case, Q0142-BX165 has a uniquely asymmetric |$\rm Ly\alpha$| halo, possibly due to source confusion. Consequently, we exclude it from most of the analysis that follows.

The KCWI data cube for Q2343-BX418 has been analyzed previously by Erb et al. (2018); here, we consider it in the context of the analysis of Q0142-BX165 above (see Fig. 21). The difference in peak SB between the major and minor axes CP2D spectra is nearly equal to that of Q0142-BX165 (⁠|$2.4 \times 10^{-18}~\mathrm{erg~s}^{-1}\, \mathrm{cm}^{-2}\, \mathrm{arcsec}^{-2}\mathring{\rm A}^{-1}$| for both). However, the spatial extent of the excess emission is significantly smaller for Q2343-BX418. The HST/WFC3, KCWI continuum, KCWI |$\rm Ly\alpha$| images, and OSIRIS-H α images all show that Q2343-BX418 comprises a single component whose centroids in the various images are consistent with one another. Despite its extreme SB asymmetry in |$\rm Ly\alpha$|⁠, Q2343-BX418 exhibits no other obviously peculiar property compared to the rest of the sample.

6 INTEGRATED LINE FLUX AND AZIMUTHAL ASYMMETRY

As shown in the previous sections, the degree of |$\rm Ly\alpha$| halo azimuthal asymmetry varies from case to case in our z ≃ 2.3 sample, but the correlation with the morphology of the central galaxy is sufficiently weak that, on average, |$\rm Ly\alpha$| haloes are remarkably symmetric and appear to be nearly independent – both kinematically and spatially – of the apparent orientation of the galaxy at the centre.

Thus far, we have treated the blueshifted and redshifted components of |$\rm Ly\alpha$| emission separately. However, the overall |$\rm Ly\alpha$| profile is expected to provide clues to the geometry and velocity field of circumgalactic H i. In this section, we compare the dependency between excess |$\rm Ly\alpha$| emission and |$W_{\lambda }(\rm Ly\alpha)$|⁠, total |$\rm Ly\alpha$| flux F_Ly α(tot), and the ratio of the total flux of blueshifted and redshifted components of |$\rm Ly\alpha$| emission [F_Ly α(blue)/F_Ly α(red)]. The integration windows used to compute the values in Table 1 are F_Ly α(tot): θ_tran ≤ 3 arcsec and −700 < Δv/(km s⁻¹) ≤ 1000; F_Ly α(blue): θ_tran ≤ 2.5 arcsec and −700 < Δv/(km s⁻¹) ≤ 0; and F_Ly α(red): θ_tran ≤ 3 arcsec and 0 < Δv/(km s⁻¹) ≤ 1000. Different integration windows were used for the two components in order to optimize the S/N of the integral; they were chosen based on the detected extent of each component in Fig. 17, in both θ_tran and Δv.

Fig. 22 examines whether or not there is a connection between major axis/minor axis asymmetry in |$\rm Ly\alpha$| flux and the overall |$\rm Ly\alpha$| halo properties mentioned above. For each pair of variables in Fig. 22, we indicate the value of the Pearson coefficient (r) and the corresponding probability p that the observed data set could be drawn from an uncorrelated parent sample. We also performed a linear regression using orthogonal distance regression in SciPy, which accounts for the estimated uncertainty in both x and y variables. As can be seen from the figure, most of the Pearson tests show no significant correlation and linear regression yields slopes consistent with zero. However, the Pearson test for the relation between F_major(blue) − F_minor(blue) and F_Ly α(tot) (middle left-hand panel of Fig. 22) yields p = 0.02, suggesting a marginally significant trend in which the asymmetry of the blueshifted component of |$\rm Ly\alpha$| favours the minor axis when |$F_{\rm Ly\alpha }(\rm tot)$| is weak, and the major axis when F_Ly α(tot) is strong.

The second relationship that stands out is F_major(red) − F_minor(red) and F_Ly α(blue)/F_Ly α(red) (bottom right-hand panel of Fig. 22), where a non-parametric test for correlation is not significant. The linear regression results in a marginally significant positive slope, indicating that as the blueshifted component of |$\rm Ly\alpha$| approaches the strength of the redshifted component, there is a tendency for excess emission along the major axis; for galaxies with F_Ly α(blue)/F_Ly α(red) ≲ 0.3 (i.e. smaller than the median value for the sample), the tendency is for excess |$\rm Ly\alpha$| emission along the minor axis.

To further explore the reliability of the correlations, we split the sample in two halves according to the overall halo properties and compared the subtracted CP2D spectra between the major and minor axes. Fig. 23 shows the result. As discussed in Section 5.5, the excess |$\rm Ly\alpha$| emission along the galaxy major axis for the red peak at Δv ≃ 300 km s⁻¹ within θ_tran ≲ 1 arcsec can be attributed to a single outlier (Q2343-BX418), which happens to fall above the median in all three quantities considered in Fig. 23 (i.e. on the right-hand panels of the figure). The subtracted CP2D spectra also indicate that essentially the entire excess |$\rm Ly\alpha$| emission for the blue peak along the minor axis – as identified earlier (Section 5.5) – is contributed by galaxies below the median |$W_{\lambda }(\rm Ly\alpha)$| and F_Ly α(tot) (top and middle left-hand panels of Fig. 23). In particular, the integrated significance within |$-700 \lt \Delta v / (\rm km~s^{-1}\,) \le -200$| and θ_tran ≤ 2|${_{.}^{\prime\prime}}$|5 exceeds 2.5σ for the |$W_{\lambda }(\rm Ly\alpha)\lt \mathrm{Median}$| bin. Comparison of the top two panels also illustrates the same trend of F_major(blue) − F_minor(blue) versus |$W_{\lambda }(\rm Ly\alpha)$| and F_Ly α correlations suggested by Fig. 22.

For the subsamples divided based on the value of F_Ly α(blue)/F_Ly α(red), the differences between the major axis and minor axis range of azimuthal angle are less significant: there is a marginally significant excess of |$\rm Ly\alpha$| emission, more noticeable in the redshifted component, where it appears to extend over the range θ_tran ≃ 1−3 arcsec, and the bin with higher F_Ly α(blue)/F_Ly α(red) appears to have a major axis excess over approximately the same range of angular distances. If Q2343-BX418 is removed from the stack of larger F_Ly α(blue)/F_Ly α(red) galaxies, the residual remains, showing that it is not attributable to a single outlier. However, none of the residuals in the bottom panels of Fig. 23 reaches a threshold of 2σ per resolution element.

In summary, a small statistical azimuthal asymmetry of |$\rm Ly\alpha$| haloes persists when the galaxy sample is divided into two according to central |$W_{\lambda }(\rm Ly\alpha)$|⁠, total F_Ly α, and the flux ratio of blueshifted and redshifted emission. Perhaps most intriguing is that galaxies with small or negative central |$W_{\lambda }(\rm Ly\alpha)$| have a tendency to exhibit excess |$\rm Ly\alpha$| emission along galaxy minor axes extending over a fairly large range of both θ_tran and velocity (⁠|$-700 \lt (\Delta v/\rm km~s^{-1}\,) \lt -200$|⁠). The apparent excess along the major axis of redshifted |$\rm Ly\alpha$| for the subsample with stronger |$\rm Ly\alpha$| emission and larger F_Ly α(blue)/F_Ly α(red), on the other hand, is confined to a smaller range of (redshifted) velocities, again roughly coincident with the typical location of the ‘red peak’.

7 FINER DIVISION OF GALAXY AZIMUTHAL ANGLES

At lower redshifts (z < 1), the covering fraction and column density of gas in various ionization stages are commonly found to be related to the orientation of the gaseous disc relative to the line of sight. Many authors have used background QSO or galaxy sightlines to detect strong Mg ii absorbers associated with galaxies at redshifts to allow the foreground galaxy orientations to be measured (e.g. Steidel et al. 2002; Bordoloi et al. 2011; Bouché et al. 2012; Schroetter et al. 2019; Lundgren et al. 2021). A common conclusion is that the sightlines giving rise to strong absorption tend to be those located at azimuthal angles corresponding to both the galaxy major and minor axes, but fewer (strong) absorbers are found at intermediate angles, 30° < ϕ < 60°. The kinematics of the absorbing gas also appear to be related to ϕ (e.g. Ho et al. 2017; Martin et al. 2019), with the broadest (and therefore strongest) systems found along the minor axis, presumably due to fast bi-conical outflows oriented perpendicular to the disc, followed by major axis sightlines sampling accreting or galactic fountain gas sharing the halo’s angular momentum and thus exhibiting disc-like rotation.

If the |$\rm Ly\alpha$| emission around z = 2–3 galaxies arises in CGM gas with properties similar to that of the low-ionization metallic absorbers at lower redshift, one might expect to see similar evidence for asymmetries along the two principal axes relative to intermediate azimuthal angles. We tested this possibility by expanding our analysis in Section 5.2 by dividing ϕ into three azimuthal bins – major axis (0° ≤ ϕ < 30°), minor axis (60° ≤ ϕ ≤ 90°), and intermediate (30° ≤ ϕ ≤ 60°). CP2D difference spectra among these three bins are shown in Fig. 24.

Fig. 24 shows no obvious sign of a bimodal distribution of |$\rm Ly\alpha$| emission with respect to ϕ, in which case the middle and right-hand panels would show residuals of opposite sign. Instead, the residual maps of ‘Major−Intermediate’ and ‘Intermediate−Minor’ suggest a gradual transition of the small asymmetries between major and minor axes in the larger bins of ϕ shown previously.

In any case, the main conclusion to draw from Fig. 24 is once again that |$\rm Ly\alpha$| emission haloes are remarkably similar in all directions with respect to the projected principal axes of z ∼ 2–3 galaxies.

8 DISCUSSION

8.1 Comparison to previous work

The morphology of |$\rm Ly\alpha$| emission from the CGM and its relation to the host galaxies has been analyzed in various works between z = 0 and z ≲ 4. However, analyses quantifying the |$\rm Ly\alpha$| emission with respect to the host galaxy orientation are limited. In this section, we attempt to compare the existing research on |$\rm Ly\alpha$| halo morphology with our findings.

At very low redshifts (z = 0.02–0.2), Guaita et al. (2015) studied the |$\rm Ly\alpha$| haloes of 14 |$\rm Ly\alpha$|-emitting galaxies in the Lyman Alpha Reference Sample (LARS; Östlin et al. 2014) and their connection to the host galaxy morphologies. They found that the |$\rm Ly\alpha$| haloes for the subset of the sample that would be considered Lyα emitters (LAE) are largely axisymmetric, and there is no single-galaxy morphological property that can be easily connected to the overall shape of the |$\rm Ly\alpha$|-emitting regions. Indeed, the |$\rm Ly\alpha$| images of the individual galaxies in Hayes et al. (2014) appear to be independent of galaxy orientation beyond D_tran ∼ 2 pkpc. Meanwhile, the stacked |$\rm Ly\alpha$| image of the LARS galaxies is elongated in the same direction as the far-UV continuum, i.e. along the major axis. However, this stack could be significantly affected by sample variance from galaxies that are particularly bright in rest-UV and |$\rm Ly\alpha$|⁠, as we saw for the KBSS sample before removing outliers. Interestingly, Duval et al. (2016) studied a special galaxy in the LARS sample, which is almost perfectly edge-on – they found two small |$\rm Ly\alpha$|-emitting components (each of extent <1 pkpc) near the disc consistent with |$\rm Ly\alpha$| emission escaping from the ISM through ‘holes’ in the galactic disc, akin to that expected in a classical Galactic fountain (Bregman 1980). While this is likely to be driven by stellar feedback allowing |$\rm Ly\alpha$| to escape in the direction perpendicular to the galaxy disc, the observed |$\rm Ly\alpha$| features are far closer to the galaxy than could be measured in our z ∼ 2–3 sample. Moreover, most of the galaxies in our sample do not have organized thin discs and so likely have much lower dust column densities obscuring active star-forming regions.

At z > 3, Bacon et al. (2017) (and subsequent papers from the same group) conducted a systematic survey of |$\rm Ly\alpha$|-emitting galaxies in the Hubble Ultra Deep Field using VLT/MUSE. For example, Leclercq et al. (2017) found significant variation of |$\rm Ly\alpha$| halo morphology among 145 galaxies and identified correlations between the halo size and the size and brightness of the galaxy UV continuum. However, its connection with the galaxy morphological orientation remains to be investigated.

Meanwhile, Chen et al. (2020a) examined a case of a strongly lensed pair of galaxies with extended |$\rm Ly\alpha$| emission at z > 3. The continuum image in the reconstructed source plane of their system A has at least three subcomponents extending over >1 arcsec, which may be similar to the subsample of galaxies shown in our Fig. 6, in terms of morphological complexity. In the context of the analysis we describe in the present work, this arrangement of |$\rm Ly\alpha$| with respect to continuum emission would be classified as excess minor axis |$\rm Ly\alpha$| emission. Unfortunately, the |$\rm Ly\alpha$| halo of a second z > 3 system is truncated in the source plane reconstruction.

In summary, we compared our result with galaxies and their |$\rm Ly\alpha$| emission morphology at z ∼ 0 and z > 3 in previous work, finding that although our results are qualitatively consistent with earlier results, the comparison is hampered by limited sample sizes, as well as by differences in redshift and intrinsic galaxy properties represented in each sample.

8.2 Theoretical predictions

Many existing studies of the distribution and kinematics of CGM gas have focused on simulated galaxies within cosmological hydrodynamic simulations. Péroux et al. (2020) analyzed how inflowing and outflowing gas might distribute differently as a function of the galactic azimuthal angle within the EAGLE and IllustrisTNG simulations, finding significant angular dependence of the flow rate and direction, as well as the CGM metallicity, with outflows of higher metallicity gas favouring the galaxy minor axis, and accretion of more metal-poor gas tending to occur along the major axis at z < 1. Although Péroux et al. (2020) focused on z ∼ 0.5 for their study, they made clear that the predicted trends would weaken significantly with increasing redshift.

At much higher redshifts (z = 5–7) around galaxies in the FIRE suite of simulations, Smith et al. (2019) found that |$\rm Ly\alpha$| escape is highly correlated with the direction of the H i outflow. Naively, one might expect that more |$\rm Ly\alpha$| would be found along the galaxy minor axis, which is the direction along which gaseous outflows would encounter the least resistance to propagation to large galactocentric radii. However, the galaxies experiencing the most active star formation at these redshifts tend to be altered on short time-scales (∼10⁷ yr) by episodic accretion, star formation, and feedback events. Thus, the direction of outflows may change on similar time-scales, while the CGM will evolve on a longer time-scale, possibly erasing any clear signatures of alignment of outflows and |$\rm Ly\alpha$| emission. For the same reason, rapidly star-forming galaxies at z ∼ 2–3, most of which have not yet established stable stellar discs, are likely to be surrounded by gas that is similarly turbulent and disordered.

Meanwhile, analytic or semi-analytic models of |$\rm Ly\alpha$| resonant scattering for idealized outflow geometries have focused primarily on the integrated |$\rm Ly\alpha$| emission profile from the entire galaxy or |$\rm Ly\alpha$| halo. Although many models account for the impact of the geometry and kinematics of gaseous outflows or accretion on the integrated |$\rm Ly\alpha$| emission profiles, there have been fewer efforts to predict the 2D spatial and spectral profiles for detailed comparison to IFU observations.

For example, Carr et al. (2018) constructed a model to predict the |$\rm Ly\alpha$| spectral morphology assuming bi-conical H i outflows with resonant scattering. The model predicts the integrated spectral profile without spatial information. In the context of the bi-conical outflow model, an integrated profile resembling our observation is predicted when the minor axis is perpendicular to the line of sight, with a large outflow having a small opening angle. However, this particular configuration would likely give rise to a highly asymmetric spatial distribution of |$\rm Ly\alpha$| emission, which we have shown is unlikely to be consistent with our observations.

Our results highlight the need for 3D models of the cool gas in the CGM around rapidly star-forming galaxies at high redshifts, prior to the development of stable disc configuration, for which the assumption of axisymmetry of outflows, at least on average, may be closer to reality. In any case, predicting the spatial and spectral properties of |$\rm Ly\alpha$| emission will require realistic treatment of the kinematics, small-scale structure, and radiative transfer of |$\rm Ly\alpha$| photons from the sites of production to their escape last scattering from the CGM towards an observer.

Gronke et al. (2016) have devised a model that assumes a two-phase CGM, composed of optically thick clumps embedded in a highly ionized diffuse ‘inter-clump’ medium. This method has been used to fit the |$\rm Ly\alpha$| profiles at multiple locations within a spatially resolved ‘|$\rm Ly\alpha$| blob’ at z ≃ 3.1 (LAB) (Li et al. 2020). More recently, Li et al. 2021, have shown that the clumpy outflow models can be applied successfully to fit multiple regions within a spatially resolved |$\rm Ly\alpha$| halo simultaneously, i.e. using a central source producing |$\rm Ly\alpha$| which then propagates through a clumpy medium with an axisymmetric outflow (see also Steidel et al. 2011, who showed that |$\rm Ly\alpha$| emission haloes similar in extent to those presented in this paper are predicted naturally, given the observed velocity fields of outflows viewed DTB to the galaxy centre and the same radial dependence of clump covering fraction inferred from absorption line studies of background objects).

8.3 Implications for H i kinematics and |$\rm Ly\alpha$| radiative transfer

A spectral profile with a dominant redshifted component of |$\rm Ly\alpha$| emission line and a weaker blueshifted component – with peaks shifted by similar |Δv| relative to the systemic redshift – is a typical signature of an expanding geometry, a central |$\rm Ly\alpha$| source function, and resonant scattering. It has also been shown that, for an ensemble of galaxies also drawn from the same KBSS redshift survey, the mean |$\rm Ly\alpha$| absorption signature measured in the spectra of background objects within D_tran ≲ 50 pkpc (θ_tran ≲ 6|${_{.}^{\prime\prime}}$|⁠) is outflow-dominated (Chen et al. 2020b). However, although the observed |$\rm Ly\alpha$| haloes and their dependence on galaxy properties may be naturally explained by central |$\rm Ly\alpha$| sources scattering through the CGM, many authors have emphasized that collisionally excited |$\rm Ly\alpha$| emission from accreting gas (i.e. gravitational cooling – see e.g. Faucher-Giguère et al. 2010; Goerdt et al. 2010; Lake et al. 2015) and in situ photoionization by the UV background and/or local sources of ionizing photons (e.g. Leclercq et al. 2020) may also contribute significantly to extended |$\rm Ly\alpha$| haloes.

Due to the complex nature of |$\rm Ly\alpha$| radiative transfer, our results cannot resolve this issue definitively. However, the fact that the stacked CP2D |$\rm Ly\alpha$| spectra show asymmetry of <2 per cent between the major and minor axes for the ensemble, combined with (1) the empirical correlation between the central |$\rm Ly\alpha$| line strength (⁠|$W_{\lambda }(\rm Ly\alpha)$|⁠) and the |$\rm Ly\alpha$| flux integrated within the entire halo, and (2) the consistently red-peak-dominated kinematics of both the central and integrated |$\rm Ly\alpha$| line, all favour scattering of |$\rm Ly\alpha$| photons produced near the galaxy centre through an outflowing, clumpy medium – at least within D_tran ≲ 30 pkpc. The remarkable statistical symmetry of the full 2D profiles, both spatially and spectrally, suggests that most of the galaxies in our sample lack persistent disc-like configurations, an inference supported also by the ubiquity of blueshifted absorption profiles in DTB spectra of similar galaxies and their lack of dependence on HST morphology (Law et al. 2012a, b). As a consequence, outflows do not behave in the manner expected for central starbursts in disc galaxies; in other words, z ∼ 2–3 galaxies on average appear to be more axisymmetric than their lower-redshift counterparts. This may have important implications for the cycling of gas and metals into and out of forming galaxies.

On the other hand, we have detected marginal (2σ) excess |$\rm Ly\alpha$| emission along the galaxy major axis for the red peak and along the galaxy minor axis for the blue peak. While most of the excess emission along the major axis can be easily explained by the relatively small sample size and the presence of one or two extreme cases, the excess blueshifted emission along the galaxy minor axis cannot be. While it is possible that the observed asymmetries indicate the prevalence of outflows along the major axis and inflows along the minor axis (i.e. the opposite of the behaviour of galaxies at z < 1), we regard such an interpretation as unlikely. An important clue may be that most of the blueshifted, minor-axis excess is contributed by galaxies with weaker than the median |$\rm Ly\alpha$| emission strength – many in that subset have central |$W_{\lambda }(\rm Ly\alpha)\lt 0$|⁠, meaning that |$\rm Ly\alpha$| photons must scatter from higher-velocity material or in directions with a more porous distribution of optically thick gas to have a high probability of escaping. Since we can observe only the photons that escape, galaxies with lower overall |$\rm Ly\alpha$| escape fractions will mean that those that do escape must take more extreme paths on average. As shown in Section 6, the galaxies with relatively weak emission also tend to have lower values of F_Ly α(blue)/F_Ly α(red), which makes smaller absolute differences in emission strength versus azimuthal angle more noticeable.

Finally, we would like to emphasize the fact that the lack of a strong statistical correlation between the morphology of extended |$\rm Ly\alpha$| emission and the galaxy orientation does not imply that the individual |$\rm Ly\alpha$| haloes are symmetric. In fact, as shown in Figs 16, 19, and 22, individual |$\rm Ly\alpha$| haloes are often asymmetric, particularly for objects with weak central |$\rm Ly\alpha$| emission. Rather, our results indicate that morphological variations of |$\rm Ly\alpha$| haloes are uncorrelated with the apparent orientation of the host galaxy starlight.

9 SUMMARY

In this paper, we have presented the first statistical results of an IFU survey of star-forming galaxies at 〈z〉 = 2.43 drawn from the KBSS and observed with the KCWI on the Keck 2 telescope. The 59 galaxies, with stellar mass and SFR typical of the full KBSS galaxy sample, comprise the subset of the KBSS–KCWI survey with both deep KCWI observations (typical exposure times of ∼5 h) and existing high spatial-resolution images from HST and/or Keck/OSIRIS. The high-resolution images were used to determine the direction of the projected major axis of the stellar continuum light of each galaxy; the KCWI IFU data cubes were used to detect spatially and spectrally resolved |$\rm Ly\alpha$| emission from the CGM around each galaxy to a limiting SB of ≲1 × 10⁻¹⁹ ergs s⁻¹ cm⁻² arcsec⁻² in the composite data (Fig. 13), enabling detection of diffuse |$\rm Ly\alpha$| emission haloes to projected distances of θ_tran ≃ 4″  (D_tran ≃ 30 pkpc).

Our major findings are summarized below:

• We introduced ‘cylindrically projected 2D spectra’ (CP2D) in order to visualize and quantify |$\rm Ly\alpha$| spectra as a function of projected galactocentric distance D_tran. The CP2D spectra are averages of spaxels over a specified range of azimuthal angle (ϕ) at a common galactocentric distance, enabling statistical analyses of |$\rm Ly\alpha$| spectral profiles and their spatial variation simultaneously. The CP2D spectra clearly show distinct redshifted and blueshifted components of |$\rm Ly\alpha$| emission that remain distinct out to projected distances of at least 25 pkpc, with rest-frame velocity extending from |$-700 \le (\Delta v_{\rm sys}/\rm km~s^{-1}\,) \lesssim 1000$| with respect to the galaxy systemic redshift, with blue and red peaks at ≃ −300 km s|$^{-1}\,$| and ≃ +300 km s|$^{-1}\,$|⁠, respectively (Section 5.2).

• We stacked the CP2D spectra of individual galaxies after aligning their continuum major axes in bins of azimuthal angle ϕ measured with respect to the major axis. By creating difference images of the CP2D projections in independent ranges of ϕ, we showed that residual differences between ‘major axis’ and ‘minor axis’ – which would reflect asymmetries in either the spatial or spectral dimension along different ranges of ϕ – are very small, with amplitude: |$\lesssim 2\times 10^{-20}~\mathrm{erg~s}^{-1}\, \mathrm{cm}^{-2}\, \mathrm{arcsec}^{-2} \mathring{\rm A}^{-1}$|⁠, corresponding to asymmetries in |$\rm Ly\alpha$| flux amounting to |$\le 2{{\ \rm per\ cent}}$| of the total, between galaxy major and minor axis directions (Section 5.3).

• We found little evidence of statistically significant asymmetry of the |$\rm Ly\alpha$| emission, except for an excess (≃ 2σ) of |$\rm Ly\alpha$| emission along galaxy major axes for the redshifted component of |$\rm Ly\alpha$| emission, with a peak near ∼+300 km s⁻¹. However, closer inspection revealed that most of the signal was caused by two galaxies with unusually asymmetric |$\rm Ly\alpha$| haloes. After discarding these outliers, another excess emission feature, with integrated significance ≃ 2σ, manifests as excess emission in the blueshifted component of |$\rm Ly\alpha$| along galaxy minor axes. This feature extends over a large range of velocity and appears to be contributed primarily by galaxies with weaker than the sample median |$\rm Ly\alpha$| emission and central |$\rm Ly\alpha$| equivalent width |$W_{\lambda }(\rm Ly\alpha)\lt 0$|⁠. The same weak-|$\rm Ly\alpha$| subsample includes many of highest M_* galaxies in the sample, as well as many of the galaxies with the smallest flux ratio between blueshifted and redshifted components (F_Ly α(blue)/F_Ly α(red)). We speculate that this asymmetry may indicate that significant azimuthal variation of |$\rm Ly\alpha$| emission morphology exists only for galaxies with the smallest |$\rm Ly\alpha$| escape fractions within the sample, i.e. weaker |$\rm Ly\alpha$|-emitting galaxies possess more developed rotational structure. Evidently, one sees only the photons managing to find rare low-N_HI holes or that scatter from the highest velocity material, both of which are more likely along the minor axis (i.e. similar to expectations based on the standard picture of bi-conical starburst-driven outflows from disc-like systems) (Sections 5.5 and 6).

• Taken together, the results show that, statistically, the |$\rm Ly\alpha$| halo around galaxies in this sample (and, by extension, the population of relatively massive star-forming galaxies at z ∼ 2–3) has remarkably little correlation – either kinematically or spatially – with the morphological distribution of stellar continuum light of the host galaxy. The observations suggest that most of the galaxies do not conform to expectations in which outflows are bi-conical and oriented along the minor axis of disc-like configurations, with accretion occurring preferentially along the major axis, suggested by observations of CGM gas in star-forming galaxies with z < 1. Instead, the lack of systematic variation in the kinematics and spatial extent with azimuthal angle of |$\rm Ly\alpha$| emission, together with the fact that F_Ly α(blue)/F_Ly α(red) is universally smaller than unity, suggests that the bulk of |$\rm Ly\alpha$| at galactocentric distances ≲30 pkpc is scattered from the inside out. The vast majority of scattered photons propagate through a scattering medium whose kinematics are dominated by outflows that are statistically symmetric with respect to the apparent morphology of a galaxy’s starlight (Section 8).

This paper marks the first attempt to understand the relationship between Ly α emission in the CGM and its host galaxy properties in the KBSS–KCWI sample. As the central DTB |$\rm Ly\alpha$| emission was shown to be correlated with the host galaxy properties for the KBSS galaxies (e.g. Trainor et al. 2015, 2016, 2019), in forthcoming work, we will utilize the CP2D spectra to further investigate the connection between the |$\rm Ly\alpha$| halo and observable properties of the host galaxies (e.g. stellar mass, SFR, SFR surface density, etc.) to understand whether galaxies at z = 2–3 significantly impact the H i distribution in the CGM and vice versa and to place additional constraints on the source functions and radiative transfer of |$\rm Ly\alpha$| in forming galaxies. With increased sample size and improved data reduction processes, we are also pushing to higher sensitivity in the stacked spectral cubes that allow us to probe the |$\rm Ly\alpha$| spectrum at larger galactocentric distances with high fidelity.

ACKNOWLEDGEMENTS

This work has included data from Keck/KCWI (Morrissey et al. 2018), Keck/OSIRIS (Larkin et al. 2006), Keck/MOSFIRE (McLean et al. 2012), Keck/LRIS-B (Steidel et al. 2004), HST/WFC3-IR, and HST/ACS. We appreciate the contribution from the staff of the W. M. Keck Observatory and the Space Telescope Science Institute.

The following software packages have been crucial to the results presented: Astropy (Astropy Collaboration 2018), the SciPy and NumPy system (Harris et al. 2020; Virtanen et al. 2020), QFitsView,⁹ CWITools (O’Sullivan & Chen 2020), Montage,¹⁰galfit (Peng et al. 2002, 2010), and DrizzlePac.¹¹

This work has been supported in part by grant AST-2009278 from the US NSF, by NASA through grant HST-GO15287.001, and the JPL/Caltech President’s and Director’s Program (YC, CS).

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. We would like to thank the anonymous referee for providing constructive feedback. We would like to acknowledge Kurt Adelberger, Milan Bogosavljević, Max Pettini, and Rachel Theios for their contribution to the KBSS survey. It is a great pleasure for us to thank Don Neill, Mateusz Matuszewski, Luca Rizzi, Donal O’Sullivan, and Sebastiano Cantalupo for their help in handling the KCWI data, and Cameron Hummels and Max Gronke for insightful discussions. YC would like to acknowledge his grandfather, Chen Yizong, who passed away during the preparation of this manuscript.

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 the National Aeronautics and Space Administration. The Observatory was made possible by the generous financial support of the W. M. Keck Foundation.

DATA AVAILABILITY

The composite CP2D spectra and the python program used to generate figures in this article are available upon reasonable request.

Footnotes

REFERENCES