The Lyman Alpha Reference Sample XIII: High-Angular Resolution 21cm HI observations of Ly$\alpha$ emitting galaxies

The Ly$\alpha$ emission line is one of the main observables of galaxies at high redshift, but its output depends strongly on the neutral gas distribution and kinematics around the star-forming regions where UV photons are produced. We present observations of Ly$\alpha$ and 21-cm HI emission at comparable scales with the goal to qualitatively investigate how the neutral interstellar medium (ISM) properties impact Ly$\alpha$ transfer in galaxies. We have observed 21-cm HI at the highest angular resolution possible (~ 3"beam) with the VLA in two local galaxies from the Lyman Alpha Reference Sample. We contrast this data with HST Ly$\alpha$ imaging and spectroscopy, and MUSE and PMAS ionized gas observations. In LARS08, high intensity Ly$\alpha$ emission is co-spatial with high column density HI where dust content is the lowest. The Ly$\alpha$ line is strongly redshifted, consistent with velocity redistribution which allows Ly$\alpha$ escape from high column density neutral medium with low dust content. In eLARS01, high intensity Ly$\alpha$ emission is located in regions of low column density HI, below the HI data sensitivity limit ($<2\times10^{20}\,$cm$^{-2}$). The perturbed ISM distribution with low column density gas in front of the Ly$\alpha$ emission region plays an important role in the escape. In both galaxies, the faint Ly$\alpha$ emission ($\sim 1\times10^{-16}$erg.s$^{-1}$cm$^{-2}$arcsec$^{-2}$) traces intermediate H$\alpha$ emission regions where HI is found, regardless of the dust content. Dust seems to modulate, but not prevent, the formation of a faint Ly$\alpha$ halo. This study suggests the existence of scaling relations between dust, H$\alpha$, HI, and Ly$\alpha$ emission in galaxies.


INTRODUCTION
The Lyα spectral line is one of the main probes of galaxies in the high-redshift Universe, playing a central role in their detection and redshift determination through either narrow-band filter or spectroscopic observations. The possibility of using this line as a tracer of the first galaxies was initially suggested by ground state, the Lyα line is a resonant transition and is absorbed and then shortly re-emitted by any neutral Hydrogen atom it encounters, resulting in a scattering of the line both spatially and spectrally. Since the restframe wavelength of the Lyα line is located in the UV range it is prone to dust absorption, which effectively removes the absorbed photons from the Lyα output of the galaxy. These interactions with the interstellar medium (ISM) of galaxies result in a complex radiative transfer of the line, that must be understood for the Lyα spectra of high-z galaxies to be interpreted correctly and to be used to their maximum potential. However, studying the details of the ISM conditions is not possible in the high-redshift Universe where only global properties can be obtained, due to the loss of resolution and depth suffered by observations. Moreover, no facilities will be built on the current horizon that will detect 21cm emission directly from galaxies at the highest redshifts.
The Lyman Alpha Reference Samples (LARS and eLARS) are two samples containing 42 local (z = 0.028 -0.18) star forming galaxies that have been assembled to study the parameters impacting the Lyα radiative transfer in galaxies Runnholm et al. 2020, Melinder et al., in prep.). Originally based on Hubble Space Telescope (HST) imaging of the Lyα, Hα and FUV continuum emission as well as HST Cosmic Origin Spectrograph (COS) spectroscopy, the number of additional observations has grown to include a tremendous amount of supporting data allowing an in-depth investigation of the galaxy properties that influence Lyα emission. These observations comprise deep optical imaging (Micheva et al. 2018), ionized gas imaging and kinematics using the PMAS (Herenz et al. 2016) and MUSE Integral Field Units (IFUs), Infrared observations using SOFIA and Herschel, molecular gas observations with the IRAM 30m telescope and the Atacama Pathfinder EXperiment (APEX) (Puschnig et al. 2020) and neutral gas imaging and kinematics with the Green Bank Telescope (GBT) and the Karl G. Jansky Very Large Array (VLA); see Pardy et al. (2014) and Le Reste et al. (in preparation).
The galaxies of the sample have high star formation rates (see Section 2) and are close in properties to Lyman Break Galaxies (LBGs), making them analogs to some high redshift galaxies . Several studies have indeed pointed to compact star forming galaxies as local Universe proxies for highredshift star-forming galaxies (e.g., Stanway & Davies 2014;Henry et al. 2015;Izotov et al. 2015Izotov et al. , 2021. The LARS and eLARS samples therefore present an unprecedented and comprehensive multi-wavelength view on the key properties driving Lyα escape that can be reasonably compared to galaxies at higher redshift. Due to the resonant nature of the line, the Lyα emission output is expected to depend strongly on a combination of the content, geometry (Giavalisco et al. 1996) and kinematics of the HI gas in the galaxy (Kunth et al. 1998;Mas-Hesse et al. 2003;Wofford et al. 2013). High HI column densities are thought to increase Lyα scattering and thus the path length of Lyα photons. In the case where dust is present, this severely increases the probability of Lyα destruction. Therefore processes disrupting the neutral ISM have been proposed to explain the escape of Lyα photons out of galaxies. Such processes include turbulence (shifting the Lyα photons out of resonance) or feedback (carving channels in the ISM that lower the column density seen by the affected Lyα photons); (Keenan et al. 2017;Bik et al. 2018;Puschnig et al. 2020). For a review on the parameters impacting Lyα escape in local galaxies, see Hayes (2019).
Neutral gas has previously been imaged at low resolution for a subset of the galaxies of the LARS sample (∼ 46 ′′ ), thus allowing the comparison of the global HI and Lyα properties of galaxies (Pardy et al. 2014). However, low-angular resolution observations of 21-cm HI trace the gas on scales that are hardly comparable with those characterising the Lyα emission at these distances. The neutral medium of the LARS galaxies was also traced with metal absorption lines (Rivera-Thorsen et al. 2015), allowing a comparison of Lyα emission and the surrounding gas in one aperture for each galaxy. UV absorption lines can be directly compared to Lyα but most lines trace both the HI and HII gas, and are not unambiguous. For example, CII and SiII are commonly used lines which have ionization potentials of 24.4 eV and 16.3 eV respectively, while HI is ionized at 13.6 eV. Moreover, absorption spectroscopy only provides information on a single line-of-sight, and thus this technique does not assess the full distribution of neutral medium throughout the galaxy, so that the exact impact of HI distribution on Lyα emission has never been observed directly. The preferred tracer of the neutral gas medium, the 21cm HI emission line, is very faint compared to e.g Hα or Lyα and so observing it at high angular resolution requires significant surface brightness and long exposure times. Additionally, resolution in the radio domain where the 21cm line is located is inherently lower than in the optical or UV bands, making multiwavelength comparisons difficult. Therefore even at the comparatively low redshifts of the galaxies in the LARS and eLARS samples, observing the HI 21cm emission remains a challenge.
This study presents the first direct comparison of the resolved HI 21-cm and Lyα emission lines in galaxies. We have obtained VLA data of the HI line at the highest angular resolution possible (∼3 ′′ beam) for two galaxies, enabling the study of the two emission lines on scales that are comparable. We also use a wealth of additional data including IFU ionized gas observations and molecular gas data that allows the tracing of star formation regions where Lyα is produced and dust absorbing Lyα. By examining the interstellar medium phases responsible for the creation, destruction and scattering of Lyα photons, this study depicts a holistic view of the Lyα radiative transfer in star-forming galaxies, and tests current theories for Lyα escape mechanisms.
We present the sample used for this study in section 2. The observations and data reduction procedures are described in section 3. In section 4 we present the HI data products and the main results obtained through the comparison of various ISM phases and Lyα emission. In section 5 the results are discussed and we provide a summary and conclusion in section 6.

SAMPLE SELECTION
The two galaxies in this study were drawn from the LARS and eLARS galaxy samples, aiming to investigate the physical properties that impact Lyα emission in local galaxies. The LARS sample contains 14 nearby starburst galaxies (0.028 ≤ z ≤ 0.18) that were selected for their strong star formation by their Hα equivalent width (EW Hα ≥ 100Å) and FUV luminosity (9.5 ≤ log(L F UV /L ⊙ ) ≤ 10.7) (seeÖstlin et al. 2014; Hayes et al. 2014, for more details). It was later extended with the eLARS sample, containing 28 additional galaxies that evenly sample the Hα equivalent width and FUV luminosity parameter spaces (Melinder et al., in prep.). Preliminary fluxes for the eLARS sample have been presented in Runnholm et al. (2020). With lower Hα equivalent width constraints (EW Hα ≥ 40Å), the eLARS galaxies are less extreme in their properties and are representative of the galaxy population dominating the local FUV luminosity function. Together the LARS and eLARS samples are suited to comprehensively study the parameters that impact the radiative transfer of Lyα photons in local galaxies.
We have obtained VLA D and C configuration observations for nearly all the galaxies in the samples and B configuration for 15 of them; these observations will be presented in future work (Le Reste et al., in prep.). While the LARS and eLARS samples were selected to contain galaxies comparable to those identified in UV-based surveys at high-redshift (e.g. Cowie et al. 2011), the targets selected for the highest angular resolution 21cm emission followup (VLA A configuration) make up a small subset of only two galaxies: LARS08 and eLARS01. Basic properties of the galaxies in the subsample are presented in Table 1. In order to obtain the highest quality HI data with long baselines, we began our campaign targeting the galaxies with the largest 21cm flux integrals based upon VLA D and C imaging. This selection also leads to targeting more massive and gas-rich galaxies, which may also be more metal and dust-rich and thus have lower Lyα output. This can be seen from the data presented in Table 1, which shows that while the luminosities and EWs are sufficient for them to be identified in high-z Lyα surveys, the global Lyα escape fractions are rather low at ∼ 1 %. Nevertheless, we stress that the emission physics, which depends significantly upon HI and dust will be identical in all galaxies.

VLA Observations
The galaxies in the sample were observed with the D, C, B and A configurations of the Very Large Array 1 (project numbers 13A-181, 14A-077, 17A-240 and 18A-095). Observations made use of the L-band centered on the galaxies respective redshifted 21-cm HI frequency. Information about the observations for each array configuration can be found in Table 2. Lower angular resolution observations of LARS08 were previously presented in Pardy et al. (2014). We designed our observations to reach the highest possible angular resolution by following the guidelines spelled out by the NRAO 2 . These guidelines specify integration times should increase by a factor of ∼3 when stepping through the longer array configurations to reach homogeneous sensitivity limits; since we are concerned with HI emission, such common limits would correspond to a fixed gas column density. Since designing the observations we have become aware of the fact that these guidelines may be most appropriate for relatively compact emission, while in our case some of the HI gas -especially the warm neutral phase, tidal debris, and circumgalatic gas -may be extended. Nevertheless, we believe this strategy is still appropriate for several reasons. Firstly, we do not attempt to use the A and B configuration data to study the large-scale emission for which we rely solely upon  imaging of the D+C data which has higher sensitivity. The observed galaxies are at distances of 129 and 168 Mpc, which are large compared to the distances of most galaxies with interferometric HI observations. Thus our physical resolution elements are comparatively coarse (2-3 kpc at half power) in comparison to the expected size of HI clouds. Thus we expect these structures to appear compact and point-like, which would make the choice of integration times appropriate. We stress that the A-configuration data is used only when we make the comparisons on the smallest scale (i.e. with UV/optical data) and that in these cases we are likely sensitive only to the highest column density, clumpiest gas. The VLA HI data was reduced with the casa software using standard prescriptions. Radio-Frequency Interference (RFI) were removed from the data to avoid the introduction of imaging artifacts in the final data cubes. A two step procedure was applied to identify and exclude RFI. First, the automatic RFI flagging task tfcrop was applied with parameters maxnpieces=3 timecutoff=3.0 and freqcutoff=3.0. This led to the removal of ∼1% of the interferences in the timefrequency plane, consisting mainly of strong, frequencynarrow RFI. This was followed by eye inspection (using the plotms and viewer tasks) and hand-masking of the measurement sets using the flagdata task. The resulting RFI-free measurement sets were then calibrated using 1331+305=3C286 as the bandpass calibrator and either J1254+1141 or J1035+5038 as the phase calibrator. The measurement sets containing all sources were split using the task split to keep only the visibilities corresponding to the target galaxies. The calibrated targetonly measurement sets were continuum subtracted by fitting a polynomial of order 1 on line-free channels of each individual dataset using the task uvcontsub. The visibilities of each dataset were reweighted according to their scatter using the task statwt on the same channels as those where uvcontsub was applied, a necessary step when combining data from different array configurations. Typically, we used 200 line-free channels on each side of the HI line, making sure not to include channels at the ends of the bandpass.

Cleaning
The resulting datasets were cleaned to 0.5σ using the CASA task tclean in interactive mode with the auto-multithresh algorithm (Kepley et al. 2020) to identify regions containing emission. This was complemented by manual inspection of the clean region selected. The use of auto-multithresh allowed for the detection of a small HI object at the edge of the LARS08 field that was missed by eye identification due to the size of the object and remote location from the center of the field. It also helped speeding the clean region selection process considerably, especially in the highangular resolution cubes. The clean images were set to have a common beam in the tclean task and a Briggs weighting robust parameter of 0.5 was used for each individual configuration, combining low and high angular resolution configurations together (D+C, D+C+B and D+C+B+A). The spectral channel width of the cubes was set to 10 km s −1 . Additionally, a uv-taper of 50kλ was used for the D+C+B+A cubes to increase the flux sensitivity while keeping a reasonably high angular resolution. Table 3 summarizes the properties of each data cube.

Emission detection
After data reduction, the emission regions were isolated from the noise. For that purpose, the cubes were masked with dilated masking, using an algorithm developed for the identification of CO emission in local galaxies (Rosolowsky & Leroy 2006;Sun et al. 2018). The algorithm first identifies regions above a certain signal-tonoise ratio using an intensity threshold criterion snr hi and a criterion on the number of consecutive channels n hi for which the intensity threshold must be seen. This first mask is then extended spatially and spectrally to include regions with a lower intensity threshold snr lo found across a certain number of consecutive channels n lo . This type of emission detection method allows for fainter emission regions to be detected, assuming that they are connected to the higher intensity emission regions, and also limits the noise included in the final data products. After conducting tests on several datasets, we found that the parameter combination (snr hi =3.5-3, n hi =3, snr lo =1.5-2, n lo =1) and extending the mask spatially by 2 pixels allowed for the best emission recovery in our datasets while limiting noise inclusion. The masked cubes were primary beam corrected to obtain correct flux densities for the subsequent data products derivation.
3.3. HI data products 3.3.1. HI flux, mass, spectra To obtain the flux and mass associated with the HI emission, the masked, primary beam corrected cubes for the VLA D array configuration only were used. This choice is motivated by the D configuration having the highest surface brightness sensitivity, thus enabling the detection of HI flux on extended scales. Furthermore the nominal D configuration beam size at 21cm (∼ 40") is well-matched to the angular sizes of the galaxies as shown in Figure 1. Although it also contains information from the D array, the D+C+B+A configuration data is more sensitive to higher surface brightness HI due to the lower sensitivity of our longer baseline observations. Thus large-scale emission is resolved out, leading to a loss of flux.
The flux S HI in a spectral cube containing a certain number N of pixels i with flux density I i can be ex-pressed as: (1) with θ min and θ maj the beam minor and major axis in degrees, ℓ pix the angular length of a pixel in degrees and ∆v the channel width. The error on the flux was derived accounting for the 3% typical flux error on the primary calibrator (Perley & Butler 2017) and the noise in the cube. Assuming that the gas is optically thin, the HI gas mass can be derived from the flux using the following expression, where D is the luminosity distance: Spectra were obtained by summing the flux in each spectral channel of the masked, primary beam corrected D configuration data. We estimate the typical r.m.s. in adjacent frequency channels devoid of 21cm flux by collapsing the mask along the spectral direction, and computing the total noise level.

Moment maps
Moment maps were derived using the masked, primary beam corrected D+C+B+A cubes with the spectral cube Python package. Both cubes have an average synthesized beam size of 3. ′′ 6. Moment-0 maps (integrated emission maps) are derived using: Moment-1 maps are intensity-weighted velocity maps, corresponding to the velocity centroid across the spectral axis: Moment-2 maps are linked to the velocity dispersion across the spectral axis, and can be expressed as: Velocity dispersion maps can be derived from the moment-2 maps using σ = √ M 2 , although this method can sometimes lead to overestimating the true velocity dispersion (Mogotsi et al. 2016).

Additional data
3.4.1. HST data HST Lyα images were obtained using the Solar Blind Channel (SBC) of the Advanced Camera for Surveys Note-Notes on columns, from left to right: (1) -Galaxy ID.
(2) -Array configuration data used to create the cube.
(3) -Synthesized beam. (4) -Cube rms. (5) -Limiting column density. These correspond either to the 2σ or 1.5σ limit, depending on the snr lo level specified during masking. (ACS), using the F125LP and F140LP filters to isolate Lyα and the adjacent continuum, respectively. For details of the observations and data-processing, see Hayes et al. (2014) andÖstlin et al. (2014); the full procedure will soon be described in Melinder (2022, in prep). In short, we use two adjacent long-pass filters in the SBC to isolate Lyα and subtract the stellar continuum, using a spatially resolved spectral modeling method relying upon 3 additional optical HST images and spectral synthesis models. We also use HST spectroscopic data from the Cosmic Origins Spectrograph (COS): the LARS 08 spectrum was obtained as part of our own program and is presented in Rivera- Thorsen et al. (2015), while the eLARS 01 spectrum was presented in Leitherer et al. (2013). In both cases Lyα was observed through the G130M grating, which provides a nominal spectral resolution of R ≈ 17000 at the wavelength of Lyα.
The eLARS 01 data were re-processed with the COS pipeline following the methods described in Rivera-Thorsen et al. (2015).

MUSE data
Optical integral field spectroscopy of LARS08 was obtained from the Multi-Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) at the Very Large Telescope Unit Telescope 4 under ESO programme 0101.B-0703. Observations were obtained on 18 May 2018, under a DIMM seeing of ≈ 0.9 arcsec at full-widthhalf-maximum. We obtained four individual exposures, each of 650 seconds and rotated by 90 degrees between integrations. An independent sky frame (120 seconds) was also taken in order to improve the background subtraction. Data were reduced by the standard MUSE pipeline Version 2.8.3 (Weilbacher et al. 2020).
The reduced cube was continuum subtracted using a 151 pixel running median filter in each spaxel. We then fit the fluxes and velocity parameters for the nebular gas using the approximation of single-component Gaussian functions: for every spatial pixel, we simultaneously fit a large number of optical lines (≃ 40). We adopt a method that recovers a single common recession velocity and velocity dispersion for all spectral features (accounting for the variation of the MUSE spectral resolution as a function of observed wavelength), and only the normalization of each line (the flux) is allowed to vary for the lines individually. In practice this allows all lines to contribute to the velocity measurements at low SNR; at high SNR the strongest lines of hydrogen and [O III] dominate in the recovery of the kinematic properties, which in turn constrains the profile shape of the weaker lines. This effectively removes degeneracies in the fitting of weaker lines, and allows their flux to be better recovered.

PMAS data
PMAS Lens-Array (Roth et al. 2005) observations for eLARS01 were obtained under overcast conditions on March 31st, 2016 (DIMM seeing ∼2.0"). The observation is part of an integral field spectroscopic campaign to study the ionised gas kinematics of the combined LARS+eLARS sample (Herenz et al., in prep.). We used the 16 ′′ × 16 ′′ field sampled at 1 ′′ × 1 ′′ for each spectral pixel. Moreover, we used the backwards-blazed R1200 grating and read out the CCD unbinned along the dispersion axis. The exposure time was 1800s. No separate sky exposure was taken, as the targeted Hα emission line was not impacted by telluric emission. We reduced the data with the p3d pipeline (Sandin et al. 2010) following the same steps as in Herenz et al. (2016), except that we did not sky-subtract and flux calibrate the data, as neither of those steps is required for the extraction of kinematic information. Measuring the FWHM of the calibration arc lamp lines we obtain R ≈ 5300 as the final resolving power (or v FWHM = 57 km s −1 ) of the PMAS data at the observed wavelength of Hα. Other Hα kinematics observations of eLARS01 can be found in Sardaneta et al. (2020).

HI Properties
Moment maps of the HI emission in the D+C+B+A cubes for LARS08 (left) and eLARS01 (right) are presented in Figure 1, along with black or white contours showing the extent of the emission in the D+C cubes. Moment maps for the D, D+C, D+C+B and D+C+B+A are shown in Appendix Figures 8 and 9. We present channel maps for the D+C+B+A data in Figures 10 and 11. The extent of the 21cm emission differs when probed by the different array configurations. This is due to the most compact configuration data (D,C; low-resolution data) having better surface brightness sensitivity, and thus being better suited at detecting faint emission. Additionally, large-scale emission is resolved out by the most extended configuration arrays (B,A) due to the lack of short spacing between antennas. The combination of different array configuration data thus yields a compromise between brightness sensitivity and angular resolution that is sufficient to image the high surface brightness 21cm emission at high angular resolution in the regions where Lyα is produced.
The HI emission of both galaxies shows signs of galaxy interaction. For LARS08, the D+C contours indicate the presence of an HI component South West of the center of the galaxy, with an extended HI feature resembling a tidal tail likely due to a recent infall. The moment maps in Figure 8 also show they are connected in velocity-space. The D+C+B+A data, which probes gas with higher column density, does not show the companion object. Inspecting optical data from SDSS (Ahumada et al. 2020) and the Pan-STARRS1 survey (Chambers et al. 2016), two extremely faint optical objects were found on the Western edge of the companion HI object, at projected angular distances of 14 arcsec and 7 arcsec from the center of the HI cloud. Due to their offset position and their faintness, we could not unambiguously identify either of these two objects as the counterpart of the HI companion. Another small, faint HI object was found in the VLA field at a distance of 13.2 arcminutes to the center of LARS08, with no detected optical counterpart. Due to the large distance of this object to the galaxy and the lack of direct evidence for its connection, we did not extend the area shown in the figure to display it, but note this suggests that LARS08 is in a group environment.
The HI emission of LARS08 traces a prominent spiral arm in the galaxy and resembles a rotating disk extending beyond the optical diameter of the galaxy. An HI hole is present in the center of the galaxy. This feature linked to the presence of star forming regions is commonly observed in galaxies, with atomic hydrogen transitioning to a molecular phase in the center of the HI disk (see e.g. Wong & Blitz 2002;Blitz & Rosolowsky 2004;Murugeshan et al. 2019). The highest HI column density region of the galaxy (N HI > 4 × 10 21 cm −2 ) is found immediately to the West of the HI hole. eLARS01 shows a merger-like HI morphology, with tails of HI extending to the North and South of the galaxy that extend on scales of ∼ 100 kpc each. The emission in the tails can be better appreciated in the D+C contours and moment maps of the lower angular resolution data presented in Figure 9. Most of the emission in the tails is not recovered with the D+C+B+A data despite the 36 hours of integration with the VLA A configuration. This is probably due to the large angular scale over which the HI emission is distributed and hence a large fraction of the HI gas is found in an extended component that is resolved out in the D+C+B+A configuration data. The decrease in flux recovered between the D configuration-only and the D+C+B+A configuration data also supports this hypothesis (see Table 3). The central part of the galaxy shows, like in the case of LARS08, a central HI hole neighboring the highest intensity region. However, careful examination reveals the presence of an absorption feature centered on the core of one of the merging galaxies (see Appendix Figure 12. Comparable HI morpho-kinematics to those of LARS08 and eLARS01 have been seen in low angularresolution observations of some of the LARS sample galaxies (see e.g. LARS03 in Pardy et al. 2014), other nearby Lyα-emitting starburst galaxies (i.e. the extended HI tails of Tol 1924-416 in Cannon et al. (2004), see also Hayes et al. (2005);Östlin et al. (2009)), and are common in local dwarf star-forming galaxies (Jaiswal & Omar 2020). This is consistent with galaxy mergers playing a role in triggering starburst episodes in galaxies. Figure 2 shows the HI spectra of the galaxies obtained from the D configuration data. The spectrum of LARS08 shows a plateau. The associated HI flux is S HI = 1.75 ± 0.09 Jy.km s −1 (with corresponding mass M HI = 1.17 ± 0.06 × 10 10 M ⊙ ), in agreement within error bars with the value of S HI = 1.8 ± 0.18 Jy.km s −1 found with the VLA D configuration presented in Pardy et al. (2014). As already noted in that study, there is a strong discrepancy between the flux values found using the VLA D configuration and the GBT which yields S HI,GBT = 3.4 ± 0.3 Jy.km s −1 . We did not find any other HI source within the area covered by the 9 arcminutes beam of the GBT, thus the difference in flux value is most likely due to the presence of low surface brightness emission that is either spatially extended Figure 2. 21cm HI emission line spectra of LARS08 and eLARS01, obtained from the masked VLA D array configuration data. The gray shaded areas show the typical error due to the flux calibration uncertainty and the cube rms. The flux in masked velocity channels was obtained from a fixed 2D aperture corresponding to the lines of sight containing emission in the cube. The vertical dotted line shows the location of the central HI velocity, the horizontal line shows the zero level. The gray area surrounding the spectrum shows the amplitude of the error due to calibrator flux uncertainty. and resolved out by the array or below the sensitivity limit of the VLA.
The spectral profile of eLARS01 is centrally peaked and we find an associated flux S HI = 4.30 ± 0.21 Jy.km s −1 (with corresponding HI mass M HI = 1.69 ± 0.08 × 10 10 M ⊙ ). A previous HI flux measurement with the single dish Nançay radiotelescope yielded S HI = 3.19 Jy.km s −1 (Martin et al. 1991). This telescope is much less sensitive than the VLA, with a typical rms of ∼1-2mJy. Error measurements on the flux were not provided in that study, however we can reasonably conclude that the lower flux value compared to the VLA is due to the large uncertainty on the flux measurement. Table 3 shows the properties of the HI cubes corresponding to the different VLA configuration combinations.
Stacked Pan-STARRS and HST images with 21-cm HI contours of LARS08 and eLARS01 are presented in Figure 3. These figures are shown for purpose of comparison between the HI and the optical emission. For both galaxies, low surface brightness features can be seen in the optical, extending beyond the HI contours. In LARS08, low surface brightness shell features in the optical indicate that some minor companion has been tidally disrupted and annexed by the main galaxy, which supports our hypothesis that the galaxy is in a group environment.

Lyα radiative transport -LARS08
In this section and the following, we analyze the multiwavelength observations to construct a picture of the Lyα radiative transfer in LARS08 and eLARS01. Here, we examine ISM phases presented in Figure 4 in turn: the Hα emitting regions, indicating where Lyα is created, the HI emitting regions, showing the transport and scattering medium of Lyα photons, and the dust, showing where Lyα is absorbed. Figure 4 a presents contours of Hα emission on the HI column density map. Hα emission traces starforming regions, and most importantly for Lyα transfer, the intrinsic production locations of Lyα radiation within the galaxy. If we were to look at a galaxy without radiative transfer effects due to HI gas or dust, the Lyα emission morphology would thus be the same as that of the Hα. We show four Hα contours, corresponding to high Hα intensity (> 6.0 × 10 −15 erg/s/cm 2 /arcsec 2 , red contours), intermediate Hα intensity (2.5 × 10 −17 erg/s/cm 2 /arcsec 2 and 200 × 10 −17 erg/s/cm 2 /arcsec 2 , orange contours) and low Hα intensity (5 × 10 −18 erg/s/cm 2 /arcsec 2 , light orange con- Figure 3. a: Pan-STARRS g,r and i stacked image (shown with asinh stretch) of LARS08 (top) and eLARS01 (bottom) with VLA D+C+B+A (cyan) and D+C (black) 21-cm HI emission contours overlaid. HI contours increase from darker to brighter. For LARS08 the contours correspond to NHI =5 × 10 19 , 1 × 10 21 and 2.5 × 10 21 cm −2 . For eLARS01 the contours correspond to NHI =1 × 10 20 , 1 × 10 21 and 3 × 10 21 cm −2 . The white box shows the extent of the images in b and c. b: HST RGB color composite image with filters F336W (blue), F438W (green), and F775W (red) for LARS08 and filters F336W (blue), F435W (green) and F814W (red) for eLARS01 c: HST color composite showing Lyα (blue), Hα (red) and FUV (green) emission (Melinder et al., in prep.) tours). Regions of high Hα intensity corresponding to the sites of strongest star formation are indicated by numbered arrows on this figure. The Hα emission down to 5×10 −18 erg/s/cm 2 /arcsec 2 outlines the spiral arm seen in HI. We have not represented all contours for the Hα emission in LARS08 to keep the Hα levels comparable between LARS08 and eLARS01, since the MUSE Hα data used for LARS08 is more sensitive than the HST data used for eLARS01. We note that the faintest MUSE Hα contours (down to 2.5×10 −18 erg/s/cm 2 /arcsec 2 ) reach the D+C HI contours for LARS08 and thus trace the full extent of the HI disk. Figure 13 in the Appendix shows, in addition to the Hα contours presented here, CO(2-1) contours tracing molecular gas and thus reservoirs of star-forming material. The peak of the CO emission is located in the HI hole, indicating that the hydrogen in this region is found in molecular form. The CO emission likely traces a gas reservoir for star formation, which is supported by the central high intensity Hα emission contours overlapping with the CO emission. Star forming region 1 coincides with the HI hole and overlaps with the CO emission present at this location. The other two strong star forming regions that correspond to high intensity Hα emission are found to the West (region 2) and South-West (region 3) of the central star forming region. These high intensity Hα emission regions are all seen to border the highest HI column density gas (N HI > 4.0 × 10 21 cm −2 ). The strongest star formation regions are linked by an intermediate intensity contour (>2.0×10 −15 erg/s/cm 2 /arcsec 2 ), that overlaps with the highest HI column density gas.

Lyα transport -HI emission
We now investigate how the Lyα emission compares in distribution and geometry to the Hα and HI emissions. Figure 4 b presents Lyα emission contours overlaid on the HI column density map, with arrows indicating regions of strongest star formation presented on Fig. 4 a for reference. They indicate where most of the Lyα emission is produced, and where it would be observed in the absence of Lyα radiative transfer effects. Three contours are overlaid on the HI map, presenting the high intensity Lyα emission (> 5 × 10 −15 erg/s/cm 2 /arcsec 2 , dark red contour), intermediate intensity Lyα emission (> 1 × 10 −15 erg/s/cm 2 /arcsec 2 , orange contour) and faint Lyα emission (> 1 × 10 −16 erg/s/cm 2 /arcsec 2 , light orange contour). These contours were selected to show the emission peaks, the intermediate Lyα emission and the faintest contours enclosed within the HST field of view.
The strongest Lyα emission (> 5 × 10 −15 erg/s/cm 2 /arcsec 2 , dark red contour) has a clumpy morphology and is found in two locations, North of star-forming region 2 and between star-forming regions 2 and 3. The strongest Lyα emission borders, but does not coincide with, the brightest HI emission regions (N HI > 4.0 × 10 21 cm −2 ). Instead, it is found in regions with HI column density N HI ∼ 3.0 × 10 21 cm −2 at the edges of HI clouds. The fact that the highest intensity Lyα emission does not overlap with the highest HI density regions could indicate that lower spatial scattering leads to stronger Lyα emission on the line of sight. High intensity emission would thus be seen in locations where Lyα photons follow a path of least resistance through the neutral gas medium. However, the HI column density values N HI ∼ 3 × 10 21 cm −2 where the emission is seen is still very high for Lyα radiation. We discuss this surprising result and potential projection effects in section 4.2.4.

Lyα destruction -dust
Figure 4 c shows the extinction E(B-V) map of LARS08, tracing the dust content of the galaxy, with Lyα contours from Fig. 4 b overlaid. Since dust absorbs Lyα photons this effectively shows the locations where Lyα photons are likely to be removed compared to where they are seen in emission. The E(B-V) is calculated using the Hα and Hβ emission lines, assuming the extinction law from Cardelli et al. (1989) and an intrinsic line ratio of 2.86, corresponding to ionized gas with an electron density of 100 cm −3 and temperature of 10,000 K. Although the extinction map does not cover the full extent of the HI disk due to low signal-to-noise of the weaker and more obscured Hβ line in the outer part of the galaxy, some notable features can be seen. First, the high and intermediate intensity Lyα emission overlaps with regions of lower extinction and dust content (E(B- V) 0.6), consistent with the picture where the presence of dust-rich gas would lead to Lyα absorption. However surprisingly, the faint Lyα emission is seen throughout the galaxy, regardless of the dust content. This suggests that regardless of extinction, sufficient spatial scattering will allow Lyα emission to be eventually observed on the line of sight. Dust thus seems to modulate, but not prevent, extended Lyα emission.

3D analysis -the view from the COS spectra
Although the morphological comparison of Lyα to the other phases of the ISM has brought insight into what impacts the spatial scattering of Lyα photons in the galaxy, it does not solve the question of how the photons can escape from an HI medium that is extremely opaque to Lyα (N HI ∼ 10 21 cm −2 in LARS08). In this section, we make use of the HST COS spectra of the galaxies and of the 3-dimensional nature of the ionized gas IFU data and HI cubes to recover the spectral information of the various phases of the ISM presented here. Figure 5 shows the normalised HI and Lyα spectra as well as the Hα emission and the SiII λ1260 absorption as a function of the radio velocity shift in the rest-frame (using the velocity of the ionized gas, determined by the optical line fitting of the MUSE spectrum at the location of the aperture) extracted using the COS aperture 3 . We note that the COS aperture is just 2.5 arcsec in diameter, which is of the order of the synthesized beam in the D+C+B+A data. Given that Lyα is produced in the HII regions probed by Hα, this figure shows the relative velocity of the material that produces and scatters Lyα, together with the spectral profile of Lyα after the scattering process. The SiII λ1260 absorption allows one to probe the neutral medium on the same scales as Lyα and at lower column densities (N HI ≥ 10 17 cm −2 ).
We observe several important features in these spectra.
First, the SiII λ1260 absorption spectrum is blueshifted with respect to center at ∼ −140km s −1 (determined by moment-1), indicating an outflow of the neutral gas at that velocity. The absorption profile reaches 0, indicating high column density HI gas on the line of sight, and essentially complete covering. This is evidence that some of the neutral gas as seen in the 21cm HI map is indeed located in front of the Lyα emis- 3 We note that spectra presented in this figure are shifted compared to those presented for LARS08 in Rivera- Thorsen et al. (2015) due to the offset position between the SDSS aperture used to recover the redshift in that previous study and that of the COS aperture. This makes the outflow seen in LARS08 less extreme compared to what had been presented previously.
sion source, and is not seen to be overlapping due to projection effects. The HI emission line is overlapping with the Hα line, but the HI line center is slightly redshifted. At the location of the COS aperture, the Lyα intensity increases as the intensity of the HI emission decreases. The whole Lyα profile is redshifted compared to the HI and Hα lines and there is negligible emission at line-centre, indicating that the line must undergo a significant velocity redistribution on the red side to escape the scattering medium. The spectral shift between the HI and Lyα lines explain how Lyα emission can be seen to overlap spatially with high column density gas in LARS08: since Lyα is redshifted it is effectively "seeing" lower column density gas through the ISM. The Lyα is seen to peak at ∼ 200 km s −1 , with a potential secondary red peak at ∼ 330 km s −1 . While the main Lyα peak is likely produced by scattering in the red wing of the Lyα line, the potential presence of the secondary peak coincides with a value of ∼twice the expansion velocity of an ionized shell produced by the outflow, and would be expected from back-scattering on the expanding shell (Verhamme et al. 2006). There is no apparent Lyα emission on the blue side of the Lyα line, coherent with the suppression of the blue peak by neutral outflowing gas with covering fraction of 1.

Lyα radiative transport -eLARS01
Similarly to LARS08, we examine the ISM phases of eLARS01 presented in Figure 6 in turn. Figure 6 a shows Hα emission contours on the HI column density map. Numbered arrows indicate strongest Hα emission regions (>2.0×10 −15 erg/s/cm 2 /arcsec 2 ). Additionally, the right panel of Figure 13 in the appendix shows CO contours in addition to the Hα contours presented on Fig. 6 a. There are two main star forming regions in this galaxy. The high Hα intensity emission region to the North of the galaxy (region 1) is co-located with the 21cm absorption feature (see Figure  12). The second region overlaps with an HI hole and with the peak of the CO emission. This suggests that part of the HI hole corresponds to a region where the hydrogen is found in a molecular form and fuels star formation. The optical image of eLARS01 presented in Figure  3 indicates that the Hα emission regions correspond to the cores of the two galaxies that are merging. Additionally, the faintest Hα emission contour (∼2.5×10 −18 erg/s/cm 2 /arcsec 2 ) extends into the gas traced by the HI maps. Although we used a similar limiting surface brightness for our Hα contours, the size difference of the fields of view of MUSE used to get the Hα information for LARS08 and HST used for eLARS01 make it harder to quantitatively compare the Hα emission of eLARS01 to its extremely extended HI emission and to the Hα emission seen in LARS08. Figure 6 b shows Lyα emission contours overlaid on the HI column density map, with regions of strongest star-formation indicated by the same arrows as presented on Fig. 6 a. In eLARS01 the strong and intermediate Lyα emission overlaps spatially in a more direct manner with the Hα emission. This is possibly due to the lower column density HI gas in front of the Lyα photon production sites (N HI < 2.2 × 10 20 cm −2 ), leading to lower spatial scattering of the line before escape. However similarly to LARS08, the faint Lyα contour is seen to trace the Hα contour at 2.5×10 −17 erg/s/cm 2 /arcsec 2 , partly following the HI emission in the D+C+B+A data. Overall, there is a very high degree of correlation between the Hα and Lyα emission.

Lyα destruction -dust
We cannot obtain the dust distribution in eLARS01 in the same manner as for LARS08, due to the PMAS observational settings. Furthermore, E(B-V) maps obtained using HST imaging have lower fidelity due to the low SNR of Hβ emission and contamination of the Hα by the NII emission. The HST Hα/Hβ map suggests that the extinction is E(B-V)∼1 magnitude on average in the Lyα emitting region, but varies between 0.2 and 2, indicating variations in the dust geometry. In this configuration, the high intensity Lyα emission can pos- sibly escape through the regions with lower dust content within the gas. Additionally, since the HI column density is lower, Lyα does not scatter as much, and is thus less likely to be absorbed by dust. Figure 7 shows the normalised HI, Lyα and PMAS Hα emission spectra, as well as the SiII λ1260 absorption spectrum as a function of the radio velocity shift in the rest-frame extracted using the COS aperture. The 21cm spectrum extraction at the COS aperture in eLARS01 yields no measurable HI. The Lyα line peaking at 160 km s −1 is redshifted compared to the Hα emission, but this time some emission is seen at line center and a blue bump is seen at -170 km s −1 . The SiII absorption profile indicates an outflowing gas shell with expansion velocity at ∼ −290 km s −1 , with the absorption extending up to -750 km s −1 . The kinematics of the Lyα in eLARS01, in particular the smaller red peak velocity compared to LARS08 indicates that Lyα photons need less scattering to escape the neutral gas medium, certainly due to the lower gas column density. The presence of the blue bump could be due to low covering fraction of the outflowing neutral gas in eLARS01. Figure 7. Left: Comparison of the normalised HI, Lyα, Hα emission and SiII λ1260 absorption spectra as a function of the rest-frame velocity shift in eLARS01. The spectra are extracted using the COS aperture, and the rest-velocity is defined by the Hα centroid at that position. The HI spectrum is obtained from the unmasked D+C+B+A VLA data cubes. The dashed vertical line indicates the position of the rest velocity as defined by the Hα centroid in the aperture, the black horizontal line indicates the zero flux level in the case of the emission spectra, and the continuum level in the case of the absorption spectrum. Right: Zoomed-in version of the Lyα HI, and Hα intensity maps with COS aperture position represented by a black circle.

DISCUSSION
For the first time we have obtained Lyα imaging and direct observations of the atomic gas at angular resolutions where spatial comparison with Lyα becomes meaningful. We support these observations with optical integral field spectroscopy of the ionized gas obtained with the MUSE and PMAS instruments. Here we discuss the results of the comparison between different phases of the ISM and Lyα emission in the observed galaxies. We recall that the Hα emission -probed both spatially and kinematically by the MUSE and PMAS data -represents the intrinsic Lyα distribution. After production by recombination events in the HII regions, the Lyα photons then scatter in the HI material mapped by the VLA 21cm imaging. This scattering process will continue until the photons either escape, in which case we might observe them in the HST imaging if they do so in the line of sight, or are absorbed by dust, which we probe for one of the galaxies using the Hα/Hβ line ratio assuming a dust screen geometry.

The role of mergers on Lyα emission from galaxies
The HI morphologies of the two galaxies presented in this study support that both are undergoing some environmental interactions. The optical image of eLARS01 clearly shows that the two galaxies of the system are currently in a major merging event that likely triggered the starburst episode. For LARS08, shell structures in the optical image indicative of such interaction were already known (Micheva et al. 2018), but the optical evidence for the infall of a galaxy is opposite to the newly detected HI tail. The connection between LARS08 and the gas cloud without optical counterpart on the Western side of the galaxy, where the gas is metal-poor ( Figure  14) could point to the accretion of metal-poor HI gas via tidal stripping, fueling the ISM of the galaxy with neutral gas. The presence of another faint HI object detected with the VLA in the vicinity of LARS08 seem to further indicate that the galaxy is located in a group environment. HI gas is very sensitive to environmental interactions, additionally more than half of the galaxies in the LARS and eLARS sample display signs of galaxy interactions or merger in the optical. This indicates that mergers might play a decisive role in fragmenting the neutral ISM and in creating the conditions that facilitate Lyα escape.

The escape of high intensity Lyα radiation from galaxies
We have found ISM conditions helping the escape of Lyα photons in both galaxies. First, we note that high star formation rates and low HI column density seem to enable the escape of high intensity Lyα with little spatial scattering. In eLARS01, we indeed observe that the Lyα emission morphology is very similar to the Hα emission, and is coincident with a location where the column density of the HI gas overlapping with Lyα and Hα sites is below the detection limit (< 2.2×10 20 cm −2 ). Our high-angular resolution observations are particularly sensitive to high-surface brightness clouds. Thus, even if not detected with the 21cm data, it is likely that some neutral gas is present in front of the Lyα regions in a more diffuse state, or in smaller, fainter clumps. The comparison of the HI spectra also shows that the Lyα in eLARS01 needs a smaller velocity offset compared to the Hα to be able to escape the galaxy. This suggests that spatial redistribution by scattering in the surrounding gas is negligible, likely due to the presence of pathways for the Lyα emission to escape directly. Furthermore, the COS Lyα spectrum shows a blue bump peaking around v∼-170 km s −1 , which should be suppressed by the outflow of neutral gas. In order for the Lyα emission to escape and have such morphology, a relatively low dust content mixed with the HII gas in the regions where we see Lyα is also needed, although we could not evaluate this precisely with the available data. However, the Hα/Hβ line ratio obtained from HST imaging suggests a very heterogeneous dust medium, with escape through dust-free channels possible thanks to the lower HI column density which limits spatial scattering along the line of sight. High star formation rates, high HI column density, and low dust contents on the line of sight also enable the escape of high intensity Lyα emission under certain gas kinematics. In LARS08, the highest intensity Lyα emission is clumpy and is found around regions of intense star formation that border the highest column density HI clouds in regions with low dust content (Figure 4). The fact that the high intensity Lyα emission does not overlap with the highest intensity gas clouds despite the intense star formation taking place on the line of sight suggests that Lyα photons follow paths of least resistance when scattering through the gas by avoiding the densest clouds. However, it is surprising that Lyα can escape such dense medium at all. Absorption in the Lyα spectrum ( Figure 5) is visible along the line-of-sight out to velocities of almost -1000 km s −1 , which implies significant absorption of the UV continuum along the line-of-sight, with column densities in excess of 10 20 cm −2 , which is optically thick to Lyα(τ ∼ 10 5 ). Furthermore, on the line of sight where the COS spectrum is obtained, the SiII absorption line reaches zero, implying a covering fraction of 1. Figure 5 shows the Lyα is strongly redshifted with respect to rest velocity defined by nebular lines in the MUSE spectrum. The Lyα emission peaks at a velocity of ∼200 km s −1 , where co-spatial Hα decreases significantly, but where the HI density is significantly lower. Furthermore the Lyα emission as seen in the COS spectrum extends to significantly redder velocities, well beyond where detectable amounts of HI are found. Therefore, Lyα photons seem to scatter in velocity space in order to "bypass" the high column density HI gas on the line of sight and escape on the red side. We conjecture that the mechanism responsible for the velocity-space scatter of the Lyα line at the COS aperture position is a combination of scattering in the red wing of the line and back-scattering of the Lyα photons on the HI material of an expanding medium. Indeed examining the neutral gas in the COS aperture through SiII line absorption, one finds gas moving out to high velocities (−140 km s −1 on average and −365 km s −1 at maximum), guaranteeing the presence of a strong outflow. This is not seen in the HI data, however SiII probes column densities ≥ 10 17 cm −2 and can be found both in HII and HI gas, and hence is not tracing exactly the same medium as the 21-cm line is. The potential presence of a secondary peak would be expected from back-scattering of the Lyα photons on the expanding shell. The velocity redistribution of Lyα photons is probably facilitated by the low dust obscuration on this line of sight, meaning that Lyα can undergo many scattering events with limited absorption, leading to high intensity emission. On the other hand, no strong or intermediate intensity Lyα emission is seen overlapping with star forming region 1 despite the high star formation intensity, probably due to the high dust obscuration (E(B-V)> 0.6), which seems to significantly dampen the intensity of the emission at this location. Thus dust content seems to impact the escape of high intensity Lyα emission out of galaxies.

Faint Lyα emission and HI structure
We now discuss our results regarding the faint Lyα emission. We find that in both galaxies, the faint Lyα emission (1-1.5 10 −16 erg/s/cm 2 /arcsec 2 ) broadly follows the Hα emission at 5−25×10 −18 erg/s/cm 2 /arcsec 2 where HI is detected with the D+C+B+A data (N HI ≥ 2 × 10 20 cm −2 ). In LARS08, the faint Lyα traces the spiral arm, regardless of dust content on the line of sight. Thus while dust seems to impede high intensity Lyα emission, it does not prevent the formation of low intensity Lyα halo. Hayes et al. (2013) found that the size of Lyα halos in LARS galaxies was strongly correlated with quantities that also scale with dust content. This suggests that dust modulates Lyα emission by either dampening the Lyα output or allowing only the escape of Lyα emission from unobscured star forming regions, resulting in a faint Lyα halo. We note that the Lyα emission contours we can access are limited by the HST field of view. Thus our data cannot answer the question of how fainter Lyα is distributed relative to the neutral gas and whether it traces the extent of the HI in the galaxies we study here.

CONCLUSIONS AND OUTLOOK
We have presented the first comparative study of Lyα and 21-cm HI emission at high angular resolution in galaxies, supported by molecular and ionized gas observations. This study of two low-z starburst galaxies shows the observational link between Lyα creation, scattering and destruction and the local structure of the interstellar medium. The main results of our study are summarized hereafter: -Both galaxies are undergoing interactions which has possibly played a role in triggering the starburst episodes they are experiencing. In eLARS01 a major merger is ongoing, while LARS08 is interacting with a companion object with no optical counterpart. The companion is either a lowsurface brightness dwarf galaxy or a tidal remnant which has been imaged here for the first time. Additional deep optical imaging reveals shell features indicative of a prior merger.
-We find that high-intensity Lyα emission (> 5 × 10 −15 erg/s/cm 2 /arcsec 2 ) is found overlapping with regions of intense star formation traced by the Hα emission. We find two ISM states facilitating the emergence of high-intensity Lyα emission. In eLARS01, with low column density HI on the line of sight, high intensity Lyα emission can escape with little velocity redistribution and the Lyα morphology is very similar to that of the Hα emission. In LARS08, we observe that high intensity Lyα emission can escape in the presence of high column density HI with low dust extinction, but the line needs significant velocity redistribution on the red side.
-Fainter Lyα emission (1 − 1.5 × 10 −16 erg/s/cm 2 /arcsec 2 ) broadly traces regions of intermediate Hα intensity where we observe HI in the D+C+B+A data. For the galaxy where we can measure dust obscuration, it is present regardless of dust content on the line of sight, implying that increased scattering will eventually lead to the formation of a faint Lyα halo. In this context dust obscuration seems to modulate, but does not impede the formation of large-scale Lyα emission. Lower intensity Lyα emission might trace HI structures on larger scales but the limited field of view of HST compared to the 21cm data does not allow us to assess the full scale of the Lyα emission.
The results we have obtained highlight the importance of the intermediate-scale (∼kpc) ISM structure to Lyα escape, in particular high intensity Lyα emission. Although the analysis has remained very exploratory, we have been able to explain Lyα escape in both galaxies using simple scenarios that suggest that scaling relations between HI, Hα, dust and Lyα could be found.
The two galaxies presented here are interesting case studies, but they represent a very small sub-sample of the LARS+eLARS galaxies, so a quantitative analysis will require a larger sample size. Obtaining further VLA A-configuration observations for the 15 galaxies already having intermediate angular scale observations (D+C+B) would be a way to extend the results presented here to a larger range of galaxies and ISM environments. However, the use of the highest angular resolution data presents challenges in terms of exposure time and signal-to-noise ratio, that make this approach impractical for a large sample of galaxies. Additionally, we identified that the combination of D, C, and B data allows to identify several of the HI substructures and represents a compromise between angular resolution and sensitivity (see Figures 8 and 9). Thus in future work, we will exploit the VLA B-configuration data al-ready obtained for 15 LARS galaxies to explore the link identified here between Lyα escape and resolved ISM properties. Some future work will also include detailed analysis of the HI emission properties and kinematics of the full LARS + eLARS galaxies at lower angular resolution (VLA D+C data), to quantify the impact of global HI properties on Lyα emission. These will be important steps to quantify the impact of neutral gas on Lyα radiative transfer at a variety of physical scales. ACKNOWLEDGMENTS A.L.R. would like to thank the referee for their comments and members of the   Figure 12 shows the absorption profile extracted from an aperture the average size of the beam in eLARS01. The velocity offset has been calculated using the Hα centroid at the aperture location. The HI gas seen in absorption is redshifted compared to the nebular gas, indicating an inflow of neutral gas towards the star-forming region. Although a Lyα spectrum is not available at this location, the presence of the neutral gas inflow suggests that we would observe a Lyα profile with a strong blue peak at the location of 21cm absorption. It is interesting to find evidence of a neutral gas inflow overlapping with the Northern star forming region in eLARS01, just 5 kpc away from the second star forming region where the COS spectrum was obtained, and where the SiII absorption profile indicates a neutral gas outflow. This suggests that the feedback processes occurring in the two star-forming regions are different.
The absorbing column density in cm −2 can be calculated using N HI = 1.822 × 10 18 T s /f F HI,abs FHI,cont ∆v, where T s is the spin temperature in K, f is the covering fraction, F HI,abs is the absorbing flux in mJy and F HI,cont the continuum flux in mJy. We find that the continuum flux at the 21cm frequency is F HI,cont = 28.9±0.9mJy. The significant uncertainty on the spin temperature T s and the covering fraction f limits our ability to calculate the column density of the absorbing gas. On the bottom right panel of Figure 12 we show the column density as a function of the spin temperature and covering fraction for T s values characteristic of the cold neutral medium (40-500K, Kanekar et al. (2011)). We find that the absorbing gas is above the effective self-shielding limit for HI N HI > 2×10 20 cm −2 (Kanekar et al. 2011), consistent with the gas being observed in front of a star-bursting location emitting UV emission. Without self-shielding the HI at this location would likely be ionized. Since the peak of the absorption feature is found just next to the location where molecular gas is observed (see Figure 13), we conclude that the column density is likely above the threshold N HI = 5 × 10 20 cm −2 where HI starts to become molecular (Savage et al. 1977), but below N HI < 1 × 10 22 cm −2 , the threshold for which hydrogen is found in a predominantly molecular form (Schaye 2001). These limits poorly constrain the range of values for the spin temperature and the covering fraction, however the similarity between the Hα and Lyα morphologies at this location indicate that Lyα photons escape with relatively little scatter through the neutral gas of the galaxy. This could be explained by the neutral gas having a low covering fraction, which is not excluded by the limits, with covering fractions as small as 3% being possible for a spin temperature T s < 55 K.

D. CO IN LARS008 AND ELARS01
In Figure 13, we show Hα and CO(2-1) contours on the HI D+C+B+A moment-0 maps of LARS08 and eLARS01 (The CO data will be presented in Puschnig et al, in prep.). Since CO(2-1) traces molecular gas, the CO contours show star-forming regions. CO contours overlap with regions of strong star formation and low HI column density in both galaxies, suggesting the drop in HI column density is linked to the presence of Hydrogen in a molecular, rather than atomic form. E. METALLICITY MAP OF LARS08 Figure 14 presents the N2 metallicity map of LARS08, derived from the [N II ] λ6584 /Hα ratio in the MUSE data cube using the Marino et al. (2013) calibration. The metallicity exhibits a radial gradient, and is lower along the main spiral arm, especially in the outer part to the South and West of the galaxy. The strongest metallicities are found in the inner part of the galaxy, to the East in particular. Figure 11. Full HI channel map for the D+C+B+A data of eLARS01, with contours of 0.26 mJy/beam in grey and 0.51 mJy/beam in black corresponding to the 1.5σ and 3σ noise levels. The aperture used to extract the spectrum is shown by a black circle. Top right panel: 21cm absorption profile extracted from the aperture shown on the left panel. Bottom right panel: Total 21cm absorbing column density as a function of spin temperature Ts and covering fraction f for temperatures characteristic of the cold neutral medium. The solid line corresponds to the observed self-shielding limit of HI (Kanekar et al. 2011), the dotted line corresponds to the threshold for the formation of molecular hydrogen (Savage et al. 1977) and the dashed line to the threshold where most of the neutral gas is molecular (Schaye 2001). Figure 13. Zoomed-in portion of the HI column density map, with MUSE Hα contours overlaid in black and CO(2-1) contours in yellow. For LARS08 (left), the Hα contour levels increase from the outer part of the galaxy, with the same levels as in Fig.  4 a displayed. For eLARS01 (right), the Hα contour levels increase from the outer part of the galaxy, with the same levels as in Fig. 6