## Abstract

We present the first results from our CAlibrating LYMan α with Hα (CALYMHA) pilot survey at the Isaac Newton Telescope. We measure Lyα emission for 488 Hα selected galaxies at z = 2.23 from High-z Emission Line Survey in the COSMOS and UDS fields with a specially designed narrow-band filter (λc = 3918 Å, Δλ = 52 Å). We find 17 dual Hα-Lyα emitters [fLyα > 5 × 10−17 erg s−1 cm−2, of which five are X-ray active galactic nuclei (AGN)]. For star-forming galaxies, we find a range of Lyα escape fractions (fesc, measured with 3 arcsec apertures) from 2 to 30 per cent. These galaxies have masses from 3 × 108 M to 1011 M and dust attenuations E(BV) = 0–0.5. Using stacking, we measure a median escape fraction of 1.6 ± 0.5 per cent (4.0 ± 1.0 per cent without correcting Hα for dust), but show that this depends on galaxy properties. The stacked fesc tends to decrease with increasing star formation rate and dust attenuation. However, at the highest masses and dust attenuations, we detect individual galaxies with fesc much higher than the typical values from stacking, indicating significant scatter in the values of fesc. Relations between fesc and UV slope are bimodal, with high fesc for either the bluest or reddest galaxies. We speculate that this bimodality and large scatter in the values of fesc is due to additional physical mechanisms such as outflows facilitating fesc for dusty/massive systems. Lyα is significantly more extended than Hα and the UV. fesc continues to increase up to at least 20 kpc (3σ, 40 kpc [2σ]) for typical star-forming galaxies and thus the aperture is the most important predictor of fesc.

## INTRODUCTION

The Lyman α (Lyα) emission line (rest-frame 1216 Å) has emerged as a powerful tool to study distant galaxies, since it is intrinsically the brightest emission line in Hii regions and redshifted into optical wavelengths at z > 2. As a result, the Lyα line has been used to spectroscopically confirm the highest redshift galaxies (Oesch et al. 2015; Zitrin et al. 2015), select samples of galaxies with narrow-band (NB) filters (e.g. Ouchi et al. 2008; Matthee et al. 2015), find extremely young galaxies (e.g. Kashikawa et al. 2012; Sobral et al. 2015b), study the interstellar, circumgalactic and intergalactic medium (e.g. Rottgering et al. 1995; Cantalupo et al. 2014; Swinbank et al. 2015) and probe the epoch of reionization (e.g. Ouchi et al. 2010; Dijkstra 2014).

However, due to the resonant nature of Lyα, it is unknown what the observed strength of the Lyα emission line actually traces. While Lyα photons are emitted as recombination radiation in Hii regions, where ionizing photons originate from star formation or AGN activity, Lyα photons can also be emitted by collisional ionization due to cooling (e.g. Rosdahl & Blaizot 2012) and shocks. Most importantly, only a small amount of neutral hydrogen is needed to get an optical depth of 1 (with column densities of ∼1014 cm−2; Hayes 2015). Therefore, Lyα photons are likely to undergo numerous scattering events before escaping a galaxy. This increases the likelihood of Lyα being absorbed by dust and also leads to a lower surface brightness (SB; detectable as Lyα haloes; e.g. Steidel et al. 2011; Momose et al. 2014; Wisotzki et al. 2015) and diffusion in wavelength space, altering line profiles (e.g. Verhamme et al. 2008; Dijkstra 2014).

In order to use Lyα to search for and study galaxies in the early Universe, it is of key importance to directly measure the fraction of intrinsically produced Lyα (the Lyα escape fraction, fesc), and to understand how that may depend on galaxy properties. Under the assumption of case B recombination radiation, fesc can be measured by comparing the Lyα flux with Hα. Hα (rest-frame 6563 Å) is not a resonant line and typically only mildly affected by dust in a well-understood way (e.g. Garn & Best 2010). Measurements of both Lyα and Hα can thus improve the understanding of what Lyα actually traces by comparing fesc with other observables as mass, dust content, kinematics, or Lyα line properties [such as the Equivalent Width (EW) and profile].

It is in principle possible to estimate the intrinsic Lyα production using other tracers of the ionizing photon production rate [i.e. other star formation rate (SFR) indicators]. However, these all come with their own uncertainties and assumptions. For example, studies using Hβ (e.g. Ciardullo et al. 2014) and UV-selected samples (e.g. Gronwall et al. 2007; Nilsson et al. 2009; Blanc et al. 2011; Cassata et al. 2015) suffer from more significant and uncertain dust corrections, and may select a population which tends to be less dusty (e.g. Oteo et al. 2015). UV-based studies are furthermore dependent on uncertainties regarding SED modelling, and on assumptions on the time-scales (UV typically traces SFR activity over a 10 times longer time-scale than nebular emission lines; e.g. Boquien, Buat & Perret 2014). Estimates using the far-infrared (Wardlow et al. 2014; Kusakabe et al. 2015) suffer from even larger assumptions on the time-scales. Remarkably though, most studies find a consistent value of fesc ∼ 30 per cent for Lyα emitters (LAEs), and lower for UV-selected galaxies, ∼3–5 per cent (e.g. Hayes et al. 2011).

Locally, it has been found that the Lyα escape fraction anti-correlates with dust attenuation (Cowie, Barger & Hu 2010; Atek et al. 2014), although the large scatter indicates that there are other regulators of Lyα escape, such as outflows (e.g. Kunth et al. 1998; Atek et al. 2008; Rivera-Thorsen et al. 2015). However, these locally studied galaxies have been selected in different ways than typical high-redshift galaxies. Green pea galaxies (selected by their strong nebular [O iii] emission) have recently been studied as local analogues for high-redshift LAEs (e.g. Henry et al. 2015; Yang et al. 2015). These studies find indications that the escape fraction correlates with H i column density, and that is also related to galactic outflows and dust attenuation. However, the sample sizes and the dynamic range are still significantly limited.

At higher redshift, it is challenging to measure the Lyα escape fraction, as Hα can only be observed up to z ∼ 2.8 from the ground, while Lyα is hard to observe at z < 2. Therefore, z ∼ 2.5 is basically the only redshift window where we can directly measure both Lyα and Hα with current instrumentation. This experiment has been performed by Hayes et al. (2010), who found a global average escape fraction of 5 ± 4 per cent. The escape fraction is obtained by comparing integrated Hα and Lyα luminosity functions (see also Hayes et al. 2011), so the results depend on assumptions on the shape of the luminosity function, integration limits, etc. Recently, Oteo et al. (2015) found that only 4.5 per cent of the Hα emitters (HAEs) covered by Nilsson et al. (2009) are detected as LAEs, indicating a similar escape fraction.

In order to increase the sample size and study dependencies on galaxy properties, we have recently completed the first phase of our CALYMHA survey: CAlibrating LYMan α with Hα. This survey combines the z = 2.23 HAEs from High-z Emission Line Survey (HiZELS) (Sobral et al. 2013) with Lyα measurements using a custom-made NB filter (see Fig. 1. The observations from our pilot survey presented here cover the full COSMOS field and a major part of the UDS field, and are described in Sobral et al. (in preparation). The aim of this paper is to measure the escape fraction for the Hα selected sources, and measure median stacked escape fractions in multiple subsets in order to understand which galaxy properties influence fesc.

Figure 1.

Filter transmission curves for the NBs used to measure Hα (NBK) and Lyα (NB392). The NB392 filter is designed to provide complete coverage of the redshifts at which HAEs can be selected in NBK, from z = 2.20 to 2.25. The Lyα emission for all our HAEs is covered even if it is shifted by ±600 km s−1. Depending on the specific redshift, the filter transmission varies between the two lines, such that Lyα is typically over-estimated with respect to Hα (see Section 3.4.2). We statistically correct for this in stacked or median measurements.

Figure 1.

Filter transmission curves for the NBs used to measure Hα (NBK) and Lyα (NB392). The NB392 filter is designed to provide complete coverage of the redshifts at which HAEs can be selected in NBK, from z = 2.20 to 2.25. The Lyα emission for all our HAEs is covered even if it is shifted by ±600 km s−1. Depending on the specific redshift, the filter transmission varies between the two lines, such that Lyα is typically over-estimated with respect to Hα (see Section 3.4.2). We statistically correct for this in stacked or median measurements.

The structure of this paper is as follows. In Section 2 we present the sample of z = 2.23 HAEs and the Lyα observations. We describe our method to measure Lyα line-flux and escape fraction and galaxy properties in Section 3, while Section 4 describes our stacking method. Section 5 presents the Lyα properties of individual galaxies. We explore correlations between fesc and galaxy properties in Section 6 and study extended Lyα emission in Section 7. Our results are compared with other studies in Section 8 and we summarize our results and present our conclusions in Section 9.

Throughout the paper, we use a Λcold dark matter cosmology with H0 = 70 km s−1 Mpc−1, ΩM = 0.3 and ΩΛ = 0.7. Magnitudes are given in the AB system and measured in 3 arcsec diameter apertures, unless noted otherwise. At z = 2.23, 1 arcsec corresponds to a physical scale of 8.2 kpc. We use a Chabrier (2003) IMF to obtain stellar masses and star formation rates.

## SAMPLE AND OBSERVATIONS

### Sample of Hα emitters

We use a sample of HAEs at z = 2.23 in the COSMOS and UDS fields selected from the HiZELS (Geach et al. 2008; Sobral et al. 2009; Best et al. 2013; Sobral et al. 2013) using NB imaging in the K band with the United Kingdom Infra Red Telescope (UKIRT). HAEs are identified using BzK and BRU colours and photometric redshifts, as described in Sobral et al. (2013). These HAEs are selected to have EW$$_{0, \rm H\alpha +[N\,\small {II}]} > 25$$ Å. In total, there are 588 HAEs at z = 2.23 in COSMOS, of which 552 are covered by our Lyα survey area. We remove 119 HAEs because they are found in noisy regions of the Lyα coverage, resulting in a sample of 433 HAEs in COSMOS. The UDS sample consists of 184 HAEs, of which 55 are observed to sufficient signal-to-noise ratio in the Isaac Newton Telescope (INT) imaging (local background of 23.5, 3σ, or deeper). This means that our total sample includes 488 HAEs, shown in Fig. 2.

Figure 2.

Positions on the sky in COSMOS and UDS of the HAEs from Sobral et al. (2013) in red points, where the size of the symbols scales with observed Hα luminosity. Our ∼2 deg2 coverage includes a wide range of environments, with number density of sources on the sky varying over orders of magnitudes, overcoming cosmic variance (see e.g. Sobral et al. 2015a). The grey points show all detections in our NB392 observations, after conservative masking of noisy regions due to the dithering pattern. It can be seen that some pointings are shallower with a lower number density of sources, and that we masked regions around bright stars and severe damages to one of the chips. After our conservative masking, we use a total area of 1.208 deg2 in COSMOS and 0.224 deg2 in UDS. We also show the four detector chips of the WFC on the INT with a total field of view of ∼0.25 deg2.

Figure 2.

Positions on the sky in COSMOS and UDS of the HAEs from Sobral et al. (2013) in red points, where the size of the symbols scales with observed Hα luminosity. Our ∼2 deg2 coverage includes a wide range of environments, with number density of sources on the sky varying over orders of magnitudes, overcoming cosmic variance (see e.g. Sobral et al. 2015a). The grey points show all detections in our NB392 observations, after conservative masking of noisy regions due to the dithering pattern. It can be seen that some pointings are shallower with a lower number density of sources, and that we masked regions around bright stars and severe damages to one of the chips. After our conservative masking, we use a total area of 1.208 deg2 in COSMOS and 0.224 deg2 in UDS. We also show the four detector chips of the WFC on the INT with a total field of view of ∼0.25 deg2.

The multi-wavelength properties of the HAEs are discussed in Oteo et al. (2015), showing that the Hα selection incorporates the full diversity of star-forming galaxies (SFGs; e.g. in their figs 5 and 6), while selections based on the Lyman break or the Lyα emission line miss significant parts of the SFG population at z = 2.23. Furthermore, although our sample of galaxies contains strongly star-bursting systems, the majority is not biased towards these rare sources. Our sample is dominated by typical galaxies which are on the main relation between stellar mass and SFR (see fig. 10 in Oteo et al. 2015, and e.g. Rodighiero et al. 2014).

### Lyα observations at z = 2.23

Lyα observations were conducted at the INT at the Observatorio Roque de los Muchachos on the island of La Palma with a specially designed NB filter for the Wide Field Camera (WFC). This NB filter (NB392, λc = 3918 Å, Δλ = 52 Å) was designed for our survey such that it observes Lyα emission for all redshifts1 at which HAEs can be selected with the NBK filter (see Fig. 1).

The details of the observations, data reduction and calibration are presented in Sobral et al. (in preparation), where we also present the Lyα luminosity function (LF), and other line-emitters detected in our NB data, such as C iv1549 at z ≈ 1.5 and [O ii] at z ≈ 0.05. For the purpose of this paper, we use the INT observations to measure the Lyα flux from Hα selected galaxies by creating thumbnail images in NB392. For continuum estimation in COSMOS, we align publicly available U and B bands (from Canada–France–Hawaii Telescope and Subaru, respectively, Capak et al. 2007; McCracken et al. 2010) and measure the flux in these filters at the positions at which the Hα emission is detected. In UDS, we use Canada–France–Hawaii Telescope U-band data (PI: Almaini & Foucaud) from UKIDSS UDS (Lawrence et al. 2007) and Subaru B-band data from SXDS (Furusawa et al. 2008).

We converted the U, B, NBK and K images to the pixel scale of the INT WFC (0.33 arcsec pixel−1). The astrometry of COSMOS images is aligned using scamp (Bertin 2006), with a reference coordinate system based on HST ACS F814W band observations (as in the public COSMOS data; McCracken et al. 2010). The UDS images are aligned to 2MASS (Skrutskie et al. 2006). The accuracy of the astrometry is of the order of 0.1 arcsec. We match the full width at half-maximum (FWHM) of the point spread function (PSF) of all images to the FWHM of the NB392 observations (ranging from 1.8 to 2.0 arcsec, depending on the particular pointing). The FWHM of reference stars was measured with SExtractor (Bertin & Arnouts 1996), which fits a Gaussian profile to the upper 80 per cent of the light profile from each detected object. For NB392 imaging, we selected reference stars with magnitudes ranging from 16 to 18, resulting in ∼20 stars per WFC detector. The reference stars in U(B) are fainter because in U(B) stars with magnitude <18(19) are saturated. For each frame, we find ∼50 reference stars with magnitudes ranging from 19 to 21. PSF matching was then done by convolving images with a Gaussian kernel. This procedure is based on the PSF matching procedure from the Subaru Suprime-Cam data reduction pipeline (Ouchi et al. 2004).

## MEASUREMENTS

### Choice of aperture

Due to resonant scattering of Lyα photons, the choice of aperture can have an important consequence on the measured Lyα flux and escape fraction, particularly given the evidence of extended Lyα emission for a range of SFGs (e.g. Steidel et al. 2011; Momose et al. 2014 and which we confirm for our sample in Section 7). Previous surveys of LAEs typically used mag-auto photometry with SExtractor to measure Lyα fluxes (e.g. Hayes et al. 2010; Ouchi et al. 2010). However, the measured flux with mag-auto will be dependent on the depth of the NB imaging. As we are measuring Lyα emission for Hα selected galaxies at the position of Hα detection, it is impossible to perform a similar mag-auto measurement as Lyα selected surveys without uncontrolled bias. This is because we have no a priori knowledge of the optimal aperture to measure Lyα. In fact, we find in Section 5 that most HAEs are undetected in Lyα at the flux limit of our observations. We also note that mag-auto measurements are dependent on the depth, and therefore are not suitable for an optimal comparison as the depth of our survey varies across the field and is different than other surveys.

Due to these considerations, we choose to use fixed diameter aperture measurements for individual sources. An aperture size of 3 arcsec was chosen for the following reasons. First, it corresponds to a radial distance of 12 kpc, which is larger than the exponential scale length of Lyα selected sources at z = 2.2–6.6 of 5–10 kpc (Momose et al. 2014), and which is also similar to the reference scale used in the study of individual Lyα haloes (Wisotzki et al. 2015; although note that this survey has detected extended Lyα emission up to a radial distance of 25 kpc). Secondly, we find that 3 arcsec aperture magnitudes on the PSF convolved images of the U, B, NBK and K band recover similar magnitudes as the 2 arcsec diameter apertures on the original Hα images (which typically have a PSF FWHM ∼0.8 arcsec), with a standard deviation of 0.2 mag. These magnitudes from 2 arcsec aperture measurements are used on most studies of the HAEs from our sample (e.g. Sobral et al. 2013; Oteo et al. 2015). For stacks of subsets of HAEs we vary the aperture, and discuss the difference in Section 7.

### Measuring line-fluxes

We use fluxes in NB392, U and B band to measure the Lyα line-flux on the positions of the HAEs using dual-mode SExtractor. The NB392 flux is calibrated on U-band magnitudes of photometrically selected galaxies (see Sobral et al., in preparation), since stars have the strong Caii3933 absorption feature at the wavelengths of the NB filter. After this calibration, we also make sure that the NB excess (U − NB392) is not a function of the UB colour, such that a very blue/red continuum does not bias line-flux measurements. This means that we empirically correct the NB magnitude using:

(1)
$${\rm NB}392_{\rm corrected} = {\rm NB}392 + 0.19\times (U-B)-0.09.$$
This correction ensures that a zero NB excess translates into a zero line-flux in NB392. For sources which are undetected in U or B we assign the median correction of the sources which are detected in U and B, which is +0.02. In the following, we refer to the broad-band U as BB. Then, with the NB and continuum measurements, the Lyα line-flux is calculated using:
(2)
$$f_{\rm Ly\alpha }= \Delta \lambda _{\rm NB} \frac{f_{\rm NB}-f_{\rm BB}}{1-\frac{\Delta \lambda _{\rm NB}}{\Delta \lambda _{\rm BB}}}.$$
Here, fNB and fBB are the flux-densities in NB392 and U and ΔλNB and ΔλBB the filter-widths, which are 52 Å and 758 Å, respectively.

We measure Hα line-fluxes as described in Sobral et al. (2013). The relevant NB is NBK and the continuum is measured in K band. The excess is corrected with the median correction of +0.03 derived from HK colours. For an HAE to be selected as a double Hα-LAE, we require the U − NB392 excess to be >0.2 (corresponding to EW0 > 4 Å) and a Lyα excess significance Σ > 2 (cf. Bunker et al. 1995; Sobral et al. 2013), see the dashed lines in Fig. 3, which we base on local measurements of the NB and broad-band background in empty 3 arcsec diameter apertures. This relatively low excess significance is only appropriate because we observe pre-selected HAEs. We note, however, that all our directly detected sources are detected with at least 3σ significance in the NB392 imaging.

Figure 3.

NB excess diagram of the sources in COSMOS and UDS. Grey points show all NB392 detections, where U has been measured in dual-mode. The green points show the HAEs which are directly detected in the NB392 imaging, with measurements done at the position of the HAEs. The red triangles are upper limits at the positions of the HAEs. The blue horizontal lines show to which rest-frame Lyα EW a certain excess corresponds. Dashed black lines show the excess significance for either the shallowest (left) or deepest (right) NB392 data. Note that some upper limits on the NB392 magnitude are actually weaker than some detections. This is due to variations in the depth of our NB392 observations across the field. Many stars have a negative excess due to the Caii3933 absorption feature.

Figure 3.

NB excess diagram of the sources in COSMOS and UDS. Grey points show all NB392 detections, where U has been measured in dual-mode. The green points show the HAEs which are directly detected in the NB392 imaging, with measurements done at the position of the HAEs. The red triangles are upper limits at the positions of the HAEs. The blue horizontal lines show to which rest-frame Lyα EW a certain excess corresponds. Dashed black lines show the excess significance for either the shallowest (left) or deepest (right) NB392 data. Note that some upper limits on the NB392 magnitude are actually weaker than some detections. This is due to variations in the depth of our NB392 observations across the field. Many stars have a negative excess due to the Caii3933 absorption feature.

### Measuring the Lyα escape fraction

In order to measure the observed fraction of Lyα flux, we need to carefully estimate the intrinsic Lyα line-flux. The intrinsic emission of Lyα due to recombination radiation is related to the Hα flux, and scales with the number of ionizing photons per second. Assuming case B recombination, a temperature T between 5000 and 20 000 K and electron density ne ranging from 102 and 104 cm−3, the intrinsic ratio of Lyα/Hα ranges from 8.1 to 11.6 (e.g. Hummer & Storey 1987). For consistency with other surveys (as discussed by e.g. Hayes 2015; Henry et al. 2015), we assume ne ≈ 350 cm−3 and T = 104K, such that the intrinsic ratio between Lyα and Hα is 8.7. Therefore, we define the Lyα escape fraction as:

(3)
$${f}_{\rm esc} = \frac{f_{\rm Ly\alpha , obs}}{8.7 f_{\rm H\alpha , corrected}}.$$
In the presence of an AGN, the assumption of case B recombination is likely invalid because e.g. collisional ionization might play a role due to shocks, leading to false estimates of the escape fraction. Among the sample of HAEs, we identify nine X-ray AGN in COSMOS using Chandra detections (Elvis et al. 2009), which are all spectroscopically confirmed to be at z = 2.23 (Civano et al. 2012). Eight of these are significantly detected in NB392 imaging. We exclude these AGN from stacking analyses, but will keep them in our sample for studying individual sources.

Note that since we measure line-fluxes in 3 arcsec apertures, fesc is strictly speaking the escape fraction within a radius of 12 kpc. It is possible that the total escape fraction is higher, particularly in the presence of an extended low SB halo due to resonant scattering (see also the discussion from a modeller's point of view by Zheng et al. 2010).

### Corrections to measurements

Although our matched NB survey requires less assumptions and uncertain conversions than escape fraction estimates based on UV or other emission-line measurements, we still need to take the following uncertainties/effects into account:

• interlopers in the Hα sample (Section 3.4.1)

• different filter transmissions (Section 3.4.2)

• dust correction of the observed Hα flux (Section 3.4.3)

• [Nii] contributing to the flux in the Hα filter (Section 3.4.4).

#### Interlopers

Our Hα sample is selected using photometric redshifts and colour–colour techniques in a sample of emission line galaxies obtained from NB imaging (see Sobral et al. 2013). This means that galaxies with other emission lines than Hα can contaminate the sample if the photometric redshift is wrongly assigned (e.g. if the galaxy has anomalous colours). Spectroscopic follow-up shows a 10 per cent interloper fraction, although this follow-up is so far limited to the brightest sources. These interlopers are either dusty low-redshift (z < 1) sources, such as Paβ at z = 0.65, or Hβ/[O iii] emitters (z ∼ 3.2–3.3; e.g. Khostovan et al. 2015). For the z ∼ 3.3 emitters, the NB392 would only measure noise, as the NB392 filter observes below the Lyman break for higher redshift galaxies and the flux for the low-redshift interlopers is typically much fainter than the NB392 limit. The identified interlopers do not occupy a particular region in the parameter space of the sample of HAEs. There may be small dependences of contamination with galaxy properties, but no trends are seen for our limited follow-up, thus we assume a flat contamination. For stacking, we increase our observed NB392 flux by 10 per cent to account for these interlopers. For individual sources without NB392 detection, we are careful in our analysis as there is the risk of interlopers, even though the fraction is relatively small.

#### Filter transmissions

While the NB392 and NBK filters are very well matched in terms of redshift coverage, the transmission at fixed redshift varies between Hα and Lyα. This means that the measured escape fraction is influenced by the particular redshift of the galaxy and resulting different filter transmissions for Hα and Lyα. Furthermore, systematic velocity offsets between Lyα and Hα might increase this effect, as it has been found that Lyα is redshifted typically 200 (400) km s−1 with respect to Hα in Lyα (UV) selected galaxies (e.g. Steidel et al. 2010; McLinden et al. 2011; Kulas et al. 2012; Hashimoto et al. 2013; Erb et al. 2014; Shibuya et al. 2014; Song et al. 2014; Sobral et al. 2015b; Trainor et al. 2015). We test the effect of the different filter transmissions and velocity shifts using a Monte Carlo simulation, similar to e.g. Nakajima et al. (2012). We simulate 1000 000 galaxies with redshifts between the limits of the NBK filter, and with a redshift probability distribution given by the NBK filter transmission (as our sample is Hα selected). Then, we redshift the Lyα line w.r.t. Hα with velocity shifts ranging from 0 to 800 km s−1, and fold it through the filter transmission in NB392. Finally, we compute the average relative Hα–Lyα transmission. For a zero velocity offset, the average transmission is 20 per cent higher for Lyα than Hα, because the NB392 filter is more top-hat like than the NBK filter. Increasing the velocity offset leads to an average lower Lyα transmission, as it is redshifted into lower transmission regions in the right wing of the filter. This effect is, however, very small, as it is constant up to a velocity shift of 400 km s−1, and decreases to 11 per cent for 800 km s−1. Because of this, we decrease the Lyα–Hα ratio of stacks and individual sources by 20 per cent. We add the 20 per cent uncertainty of this correction to the error on the escape fractions in quadrature. Spectroscopic follow-up is required to fully investigate the effect of velocity offsets on our measured escape fractions.

Table 1.

Numbers and median properties with 1σ deviations of the sample of Hα-LAE with and without AGN. We also show the upper limits on the galaxies that are not detected in Lyα, which is the comparison sample. Note that these are the median upper limits. Hα and Lyα fluxes are the values observed in 3 arcsec apertures. For completeness, we show the subsample of SFGs that we used for stacks (these are selected based on the depth of Lyα observations). The masses are derived from SED fitting (Sobral et al. 2014), which also gives Hα attenuation based on the stellar extinction and the Calzetti et al. (2000) law (see Section 3.5). The Hα attenuation from Garn & Best (2010) is based on a calibration between dust, SFR and mass. The total sample consists of 488 HAEs, with a stacked median escape fraction of 0.9 ± 0.1 per cent (for 3 arcsec diameter apertures), which is lower than the median escape fraction of individually detected source, because for individual sources we are observationally biased towards high escape fractions. *This escape fraction is likely wrong, as in AGN there is likely a departure from case B recombination due to shocks. We still show this for comparison, indicating that Lyα is typically bright for AGN.

Hα sample Nr. f〉 fLyα〉 log10(MstarAHα,Calzetti $$A_{\rm H\alpha , Garn\&Best}$$ fesc
[10−16 erg s−1 cm−2[10−16 erg s−1 cm−2[M[mag] [mag] [ per cent]
SFG with Lyα 12 0.5 ± 0.3 0.7 ± 0.5 10.3 ± 0.8 0.83 ± 0.4 1.11 ± 0.4 10.8 ± 1.3
SFG no Lyα 468 0.4 ± 0.3 <0.7 9.9 ± 0.7 0.83 ± 0.5 0.86 ± 0.4 <20.1
AGN with Lyα 1.3 ± 0.4 3.6 ± 1.5 10.8 ± 0.4 0.50 ± 0.3 1.55 ± 0.3 12.8 ± 1.4*
SFGs for stacks 265 0.4 ± 0.3 0.1 ± 0.01 9.9 ± 0.7 1.0 ± 0.4 0.86 ± 0.4 0.9 ± 0.1
Hα sample Nr. f〉 fLyα〉 log10(MstarAHα,Calzetti $$A_{\rm H\alpha , Garn\&Best}$$ fesc
[10−16 erg s−1 cm−2[10−16 erg s−1 cm−2[M[mag] [mag] [ per cent]
SFG with Lyα 12 0.5 ± 0.3 0.7 ± 0.5 10.3 ± 0.8 0.83 ± 0.4 1.11 ± 0.4 10.8 ± 1.3
SFG no Lyα 468 0.4 ± 0.3 <0.7 9.9 ± 0.7 0.83 ± 0.5 0.86 ± 0.4 <20.1
AGN with Lyα 1.3 ± 0.4 3.6 ± 1.5 10.8 ± 0.4 0.50 ± 0.3 1.55 ± 0.3 12.8 ± 1.4*
SFGs for stacks 265 0.4 ± 0.3 0.1 ± 0.01 9.9 ± 0.7 1.0 ± 0.4 0.86 ± 0.4 0.9 ± 0.1

#### Dust attenuation

Even though the Hα emission line is at red wavelengths compared to, for example, UV radiation, it is still affected by dust, such that we underestimate the intrinsic Hα luminosity. Correcting for dust typically involves a number of uncertainties, such as the shape and normalization of the attenuation curve, the difference between nebular and stellar extinction (e.g. Reddy et al. 2015 and references therein) and the general uncertainties in SED fitting. For consistency with other surveys, we correct for extinction by applying a Calzetti et al. (2000) dust correction, using the estimated extinction, E(BV)star, measurements from the best-fitting SED model from Sobral et al. (2014). Note that we assume E(BV)star = E(BV)gas, independent of galaxy property. Recent spectroscopic results at z ∼ 2 (e.g. Reddy et al. 2015) indicate that this is reasonable when averaged over the galaxy population, although there are indications that the nebular attenuation is higher than the stellar attenuation for galaxies with high SFR, particularly for galaxies with SFR > 50 M yr−1. Therefore, if such a trend would be confirmed, our inferred relations between fesc may be slightly affected. We discuss this when relevant in Sections 6 and 8.2.1.

When stacking, we use the median dust correction of the sources included in the stacked sample, which is A = 1.0. This number has also been used, for example, by Sobral et al. (2013) in order to derive the cosmic star formation rate density, which agrees very well with independent measures. Ibar et al. (2013) showed this median attenuation also holds for a similar sample of HAEs at z = 1.47 by using Herschel data.

However, we also investigate how our results change when using the dust correction prescription from Garn & Best (2010), which is a calibration between dust extinction and stellar mass based on a large sample of spectroscopically measured Balmer decrements in the local Universe. This relation between Balmer decrement and stellar mass is shown to hold up to at least z ∼ 1.5 (e.g. Sobral et al. 2012; Ibar et al. 2013).

For individual sources, the two different dust corrections explored here can vary by up to a factor of 5, as seen in Table 2. This results in large systemic errors which can only be addressed with follow-up spectroscopy to measure Balmer decrements. Throughout the paper, we add the error on the dust correction due to the error in SED fitting in quadrature to the error of the Hα flux, but note that the systematic errors in the dust correction are typically of a factor of 2.

#### [N ii] contamination

Due to the broadness of the NBK filter used to measure Hα, the adjacent [N ii] emission line doublet contributes to the observed line-flux. We correct for this contribution using the relation from Sobral et al. (2012), who calibrated a relation between [N ii]/([N ii]+Hα) and EW$$_{0,\rm H\alpha +[N\,\small {II}]}$$ on Sloan Digital Sky Survey galaxies. More recently, Sobral et al. (2015a) found the relation to hold at least up to z ∼ 1. At z = 2.23, we use this relation to infer a typical fraction of [N ii]/([N ii]+Hα) = 0.17 ± 0.08, which is consistent with spectroscopic follow-up at z ∼ 2 (Swinbank et al. 2012; Sanders et al. 2015). We have checked that our observed trends between fesc and galaxy properties do not qualitatively depend on this correction – if we apply the median correction to all sources, the results are the same within the error bars. We add 10 per cent of the correction to the error in quadrature. For stacks, we measure the EW$$_{0, \rm H\alpha +[N\,\small {II}]}$$ and apply the corresponding correction, which is consistent with the median correction mentioned here, and we also add 10 per cent of the correction to the error in quadrature.

### Definitions of galaxy properties

We compare fesc with a range of galaxy properties, defined here. SFRs are computed from Hα luminosity, assuming a luminosity distance of 17746 Mpc (corresponding to z = 2.23 with our cosmological parameters) and the conversion using a Chabrier (2003) IMF:

(4)
$${\rm SFR(H}\alpha {\rm )/(M}_{{\odot }}{\rm yr^{-1})} = 4.4\times 10^{-42}\, {L({\rm H}}\alpha {\rm )/(erg\, s}^{-1}{\rm )}$$
where L(Hα) is the dust-corrected Hα luminosity and SFR(Hα) the SFR.

Stellar masses and extinctions (E(BV)) are obtained through SED fitting as described in Sobral et al. (2014). In short, the far-UV to mid-infrared photometry is fitted with Bruzual & Charlot (2003) based SED templates, a Chabrier (2003) IMF, exponentially declining star formation histories, dust attenuation as described by Calzetti et al. (2000) and a metallicity ranging from Z = 0.0001 to 0.05. While we use a mass defined as the median mass of all fitted models within 1σ of the best fit, we use the E(BV) value of the best-fitted model. The errors on stellar masses and extinctions are computed as the 1σ variation in the fitted values from SED models that have a χ2 within 1σ of the best-fitted model. For stellar mass, these errors range from 0.2 dex for the lowest masses to 0.1 dex for the highest masses. The typical uncertainty on the extinction ranges from 0.12 at E(BV) ≈ 0.1 to 0.05 at E(BV) ≈ 0.3.

The UV slope β (which is a tracer for dust content, stellar populations and escape of continuum ionizing photons; e.g. Dunlop et al. 2012) is calculated using photometry from the observed g+R colours. These bands were chosen such that there is minimal contribution from Lyα to the g+ band (the transmission at the corresponding wavelength is <5 per cent), and such that we measure the slope at a rest-frame UV wavelength of ∼1500 Å. We also chose to derive the slope from observed colours in stead of using the SED fit, as otherwise there might be biases (e.g. the SED-based extinction correction is related to the UV slope). The error in β due to measurement errors in g+ and R ranges from typically 0.5 at β = −2.3 to 0.3 at β > 0.

## STACKING METHOD

In order to reach deeper Lyα line-fluxes, we use stacking methods to combine observations of our full sample of observed galaxies, such that the exposure time is effectively increased by a factor of ∼400. This however involves some complications and assumptions. For example, we will use the median stacked value, rather than the mean stacked value, such that our results are not biased towards bright outliers. However, our results will still be biased towards the most numerous kind of sources in our sample. Stacking also assumes that all sources are part of a single population with similar properties – which may not always be the case, as indicated by the results in the previous section.

We divide our sample in subsets of various physical properties in Section 6 and study how these stacks compare with the results from individual galaxies. We discuss the effect of varying apertures in Section 7. The errors of the measured fluxes and resulting escape fractions in stacks are estimated using the jackknife method. The errors due to differences in the PSF of the NB and broad-band are added as a function of aperture radius (see Section 4.1). We add all other sources of systematic error (see Section 3.4) in quadrature.

We obtain stacked measurements by median combining the counts in 1 × 1 arcmin2 thumbnails in U, B, NB392, NBK and K bands of the HAEs covered in our INT observations (see Fig. 2). From the stacked thumbnails (as, for example, shown in Fig. 11), we measure photometry in various apertures at the central position (defined by the position of the NBK detection. Note that our typical astrometric errors are of the order of ∼0.1 arcsec, corresponding to ∼1 kpc). With this photometry, we obtain line-fluxes for both Lyα and Hα. The Lyα flux is corrected using UB colours, and we account for the [N ii] contribution to the NBK flux using the relation with EW from Sobral et al. (2012) (see Section 3.4.4). We also add the error due to differences in the PSF of U and NB392 to the error of the Lyα flux (see Section 4.1). We apply the median dust correction of the HAEs, which is roughly similar for using the Calzetti or Garn & Best method: A = 1.0 or A = 0.86, respectively. For our full sample of 488 HAEs, we observe a median stacked Lyα line-flux of 3.5 ± 0.3 × 10−18 erg s−1 cm−2, and an escape fraction of 0.3 ± 0.06 per cent in 3 arcsec apertures, corresponding to a radial distance to the centre of ∼20 kpc. The significance of these detections is discussed in Section 7.

The depth of our NB392 observations is inhomogeneous over the full fields (see Fig. 2). We therefore study the effect of limiting our sample based on the depth of the NB392 observations. We find that the photometric errors on the stacked NB392 image are minimized when we only include sources for which the local 3σ depth is at least 24.1 AB magnitude, which corresponds to the inclusion of 265 out of the 488 sources. For the remainder of this section, we only include sources which are among these 265. The median SFRs, stellar masses, and dust attenuations of this sample are similar to the average properties of the full sample (see Table 1). The 3σ depth of the NB392 stack of these 265 sources is 27.2 AB magnitude. In the case of a pure line and no continuum contributing to the NB392 flux, this corresponds to a limiting line-flux of ∼5 × 10−18 erg s−1 cm−2.

### Empirical evaluation of different PSF shapes

The NB and broad-band observations are taken with different telescopes, cameras and at different observing sites and under different conditions. Therefore, even though we match the PSF FWHM of all images, the actual shape of the PSF might vary between NB and broad-band. This might artificially influence the SB profile of line-emission estimated from the difference between the two bands. This becomes particularly important when we study stacked images of over 300 sources, where errors on the per cent level might dominate the measured signal.

We empirically evaluate the differences in NB and broad-band PSF by performing the following sanity check: we first select line-emitters in NBK imaging, which: (i) are not selected as HAEs, (ii) are not selected as higher redshift line-emitters, or (iii) do not have a photometric redshift >1 (from Ilbert et al. 2009). With this sample, we ensure that the NB392 photometry should measure (relatively flat) continuum by removing a handful of sources with an emission line in NB392. This leaves us with 245 sources, which have a similar NBK magnitude distribution as the HAE sample. The U, B and NB392 images are stacked in the exact same way as we treat the HAEs. This is used to measure the resulting line-flux and SB profile of the stack in the NB392 band (we estimate the continuum from U and B.). Although the NB392 was photometrically calibrated to the U band in 3 arcsec diameter apertures, we detect a small residual signal with a typical SB profile of central absorption (with SB ∼−2 × 10−19 erg s−1 cm−2 arcsec−2; see Fig. 4), and peaking at a radial distance of 2.5 arcsec (with a SB of ∼4 × 10−19 erg s−1 cm−2 arcsec−2). We note that at the radial aperture of 1.5 arcsec, which we used for our calibration, the integrated flux signal is consistent with zero (see Fig. 4). Corrections therefore only need to be applied for other aperture radii and SB profiles. For individually detected Lyα sources, the residual signal is at the 1–10 per cent flux level, but for stacks, it can be more important. The origin of this residual signal is likely because of differences in the inner part of the PSF, similarly as those reported by e.g. Momose et al. (2014). The uncertainty in our astrometry is of the order of 0.1 arcsec and therefore likely less important.

Figure 4.

SB profile (blue) and integrated flux (red) of the stack of the reference sample, which should have a zero line-flux at a 1.5 arcsec aperture radius by construction. We indicate the 1.5 arcsec aperture radius with a dashed black line. The inset figure shows the 2D image obtained by subtracting the BB from the NB. We find a small residual signal which has central absorption and peak at a radial distance of 2.5 arcsec. We attribute this signal due to differences in the PSF of the NB and the broad-band. In the remainder of this paper, we subtract this signal from the SB of individual sources and from stacks and add this subtraction to the error of the total flux in quadrature. The typical signal measured for individually detected HAEs is typically 10–100 times higher than the signal due to the PSF differences, but it is of the same order of magnitude as stacked measurements.

Figure 4.

SB profile (blue) and integrated flux (red) of the stack of the reference sample, which should have a zero line-flux at a 1.5 arcsec aperture radius by construction. We indicate the 1.5 arcsec aperture radius with a dashed black line. The inset figure shows the 2D image obtained by subtracting the BB from the NB. We find a small residual signal which has central absorption and peak at a radial distance of 2.5 arcsec. We attribute this signal due to differences in the PSF of the NB and the broad-band. In the remainder of this paper, we subtract this signal from the SB of individual sources and from stacks and add this subtraction to the error of the total flux in quadrature. The typical signal measured for individually detected HAEs is typically 10–100 times higher than the signal due to the PSF differences, but it is of the same order of magnitude as stacked measurements.

We thus conclude that the differences in PSF shapes of broad-band and NB have a small effect on stacked measurements, but we still take it into account by correcting all SB profiles and any aperture measurements at values other than 3 arcsec. We add the residual flux to the error of the total flux in quadrature.

Table 2.

Properties of HAEs at z = 2.23 detected at 3σ in NB392 imaging and having a Lyα emission line. The IDs correspond to the last digits of the full HiZELS IDs, which can be retrieved by placing ‘HiZELS-COSMOS-NBK-DTC-S12B-’ in front of them. The coordinates are measured at the peak of NBK (rest-frame R+Hα) emission. The observed line-fluxes, EW, dust extinction and escape fraction are measured with 3 arcsec apertures. The escape fraction (fesc) is computed under the assumptions explained in Section 3.4, thus using AHα,C. 1: SED: dust correction using SED fitted E(BV) and a Calzetti et al. (2000) law. 2: GB: Garn & Best (2010) dust correction based on stellar mass. 3: Codes in this column correspond to: 1: X-Ray AGN, 2: [O ii] line detected in NBJ, 3: [O iii] line detected in NBH.

ID R.A. Dec. Mstar fLyα EW0, Lyα f AHα,C AHα,GB fesc Note3
(J2000) (J2000) log10(Merg s−1 cm−2 Å erg s−1 cm−2 SED1 GB2 per cent
× 10−16  × 10−16
1057 10:00:39.6 +02:02:41.2 10.6$$^{+0.1}_{-0.1}$$ 0.46 82 0.42 0.17 1.43 10.8 ± 1.4
1073 10:00:44.2 +02:02:06.9 11.1$$^{+0.1}_{-0.1}$$ 0.88 55 1.24 0.50 1.77 5.1 ± 0.6
1139 10:00:55.4 +01:59:55.4 10.8$$^{+0.1}_{-0.1}$$ 4.55 63 1.03 0.17 1.56 43.7 ± 5.1 1, 2, 3
1993 10:02:08.7 +02:21:19.9 8.6$$^{+0.1}_{-0.1}$$ 1.14 68 0.26 1.49 0.29 12.8 ± 1.5
2600 10:00:07.6 +02:00:13.2 8.7$$^{+0.2}_{-0.2}$$ 1.32 80 0.53 1.00 0.29 11.4 ± 1.3
2741 10:01:57.9 +01:54:36.9 10.1$$^{+0.1}_{-0.1}$$ 3.55 14 2.24 0.67 0.99 9.9 ± 1.0
4032 10:00:51.1 +02:41:16.9 11.0$$^{+0.1}_{-0.1}$$ 0.46 28 0.34 0.83 1.67 7.3 ± 1.1
4427 10:01:19.4 +02:07:32.6 8.7$$^{+0.2}_{-0.4}$$ 0.57 12 0.58 1.16 0.29 3.8 ± 0.6
4459 10:01:43.3 +02:11:15.7 10.3$$^{+0.1}_{-0.1}$$ 1.41 93 0.41 0.50 1.18 24.8 ± 3.1
4861 10:00:03.3 +02:11:04.4 9.0$$^{+0.2}_{-0.2}$$ 0.31 16 0.28 0.0 0.34 12.4 ± 2.3
5583 10:01:59.6 +02:39:32.7 10.8$$^{+0.2}_{-0.1}$$ 1.73 244 0.41 0.50 1.54 30.4 ± 3.8
5847 10:01:12.2 +02:53:25.9 10.3$$^{+0.1}_{-0.1}$$ 0.70 60 1.03 0.50 1.11 4.9 ± 0.5
7232 10:01:05.4 +01:46:11.6 10.3$$^{+0.1}_{-0.1}$$ 0.42 15 0.74 1.33 1.11 1.9 ± 0.3
7693 09:59:49.6 +01:50:24.7 9.6$$^{+0.1}_{-0.1}$$ 0.75 10 0.84 0.67 0.65 5.5 ± 0.8
7801 10:02:08.6 +01:45:53.6 10.4$$^{+0.1}_{-0.1}$$ 2.38 1.35 0.50 1.24 12.8 ± 1.4
9274 10:00:26.7 +01:58:23.0 11.0$$^{+0.1}_{-0.1}$$ 5.06 142 1.54 0.13 1.72 32.4 ± 3.6 1, 2
9630 10:02:31.3 +01:58:16.5 9.7$$^{+0.2}_{-0.2}$$ 0.61 29 0.44 0.0 0.70 16.1 ± 2.2
ID R.A. Dec. Mstar fLyα EW0, Lyα f AHα,C AHα,GB fesc Note3
(J2000) (J2000) log10(Merg s−1 cm−2 Å erg s−1 cm−2 SED1 GB2 per cent
× 10−16  × 10−16
1057 10:00:39.6 +02:02:41.2 10.6$$^{+0.1}_{-0.1}$$ 0.46 82 0.42 0.17 1.43 10.8 ± 1.4
1073 10:00:44.2 +02:02:06.9 11.1$$^{+0.1}_{-0.1}$$ 0.88 55 1.24 0.50 1.77 5.1 ± 0.6
1139 10:00:55.4 +01:59:55.4 10.8$$^{+0.1}_{-0.1}$$ 4.55 63 1.03 0.17 1.56 43.7 ± 5.1 1, 2, 3
1993 10:02:08.7 +02:21:19.9 8.6$$^{+0.1}_{-0.1}$$ 1.14 68 0.26 1.49 0.29 12.8 ± 1.5
2600 10:00:07.6 +02:00:13.2 8.7$$^{+0.2}_{-0.2}$$ 1.32 80 0.53 1.00 0.29 11.4 ± 1.3
2741 10:01:57.9 +01:54:36.9 10.1$$^{+0.1}_{-0.1}$$ 3.55 14 2.24 0.67 0.99 9.9 ± 1.0
4032 10:00:51.1 +02:41:16.9 11.0$$^{+0.1}_{-0.1}$$ 0.46 28 0.34 0.83 1.67 7.3 ± 1.1
4427 10:01:19.4 +02:07:32.6 8.7$$^{+0.2}_{-0.4}$$ 0.57 12 0.58 1.16 0.29 3.8 ± 0.6
4459 10:01:43.3 +02:11:15.7 10.3$$^{+0.1}_{-0.1}$$ 1.41 93 0.41 0.50 1.18 24.8 ± 3.1
4861 10:00:03.3 +02:11:04.4 9.0$$^{+0.2}_{-0.2}$$ 0.31 16 0.28 0.0 0.34 12.4 ± 2.3
5583 10:01:59.6 +02:39:32.7 10.8$$^{+0.2}_{-0.1}$$ 1.73 244 0.41 0.50 1.54 30.4 ± 3.8
5847 10:01:12.2 +02:53:25.9 10.3$$^{+0.1}_{-0.1}$$ 0.70 60 1.03 0.50 1.11 4.9 ± 0.5
7232 10:01:05.4 +01:46:11.6 10.3$$^{+0.1}_{-0.1}$$ 0.42 15 0.74 1.33 1.11 1.9 ± 0.3
7693 09:59:49.6 +01:50:24.7 9.6$$^{+0.1}_{-0.1}$$ 0.75 10 0.84 0.67 0.65 5.5 ± 0.8
7801 10:02:08.6 +01:45:53.6 10.4$$^{+0.1}_{-0.1}$$ 2.38 1.35 0.50 1.24 12.8 ± 1.4
9274 10:00:26.7 +01:58:23.0 11.0$$^{+0.1}_{-0.1}$$ 5.06 142 1.54 0.13 1.72 32.4 ± 3.6 1, 2
9630 10:02:31.3 +01:58:16.5 9.7$$^{+0.2}_{-0.2}$$ 0.61 29 0.44 0.0 0.70 16.1 ± 2.2

## DIRECT MEASUREMENTS FOR INDIVIDUAL GALAXIES

We directly detect (>3σ) 43 out of our 488 HAEs in the NB392 imaging, which is a combination of UV continuum and Lyα line (see Table 1). The 3σ limit corresponds to limiting Lyα fluxes ranging from 3.8 to 7.4 × 10−17 erg s−1 cm−2 (assuming the typical continuum level of 0.23 μJy, ∼25.5 AB magnitude in the U band). Out of these robust detections, 17 show a significant Lyα line detection (excess significance Σ > 2), all in COSMOS. The properties of these sources and their IDs from the HiZELS catalogue (Sobral et al. 2013) are listed in Table 2. The other 26 robust NB392 detections are HAEs with strong upper limits on their Lyα flux, as we have detected the UV continuum in the NB392 filter.2

Five of the dual emitters are matched (within 3 arcsec) with an X-ray detection from Chandra. From spectroscopy with IMACS and from zCOSMOS (Lilly et al. 2009), these are all classed as BL-AGN. These AGN are among the brightest and most massive HAEs (Fig. 5): all have stellar masses above 1010.5 M (Fig. 6), and the fraction of BL-AGN is consistent with the results from Sobral et al. (2016). The ISM conditions surrounding the AGN might lead to other ionizing mechanisms than case B photo-ionization, such as shocks, making it more challenging to measure the Lyα escape fraction as the intrinsic Hα-Lyα changes.

Figure 5.

Histogram of Hα fluxes in our galaxy sample at z = 2.23. AGN are typically found among the brightest HAEs, and are also typically detected in NB392, such that they are either bright in the UV continuum or Lyα or both. We also show the distribution of Lyα detected SFGs. It can be seen that not all brightest HAEs have been detected, indicating very low escape fractions or interlopers. On the other hand, some very faint HAEs are still detected in Lyα, indicative of high Lyα escape fractions.

Figure 5.

Histogram of Hα fluxes in our galaxy sample at z = 2.23. AGN are typically found among the brightest HAEs, and are also typically detected in NB392, such that they are either bright in the UV continuum or Lyα or both. We also show the distribution of Lyα detected SFGs. It can be seen that not all brightest HAEs have been detected, indicating very low escape fractions or interlopers. On the other hand, some very faint HAEs are still detected in Lyα, indicative of high Lyα escape fractions.

Figure 6.

SFR(Hα) versus stellar mass for the observed HAEs. We obtained the SFR from dust-corrected Hα and stellar mass from SED fitting (see Section 3.5). We show the position of sources with and without Lyα, AGN with Lyα and AGN without Lyα. There is no obvious difference in the SFRs or stellar masses between sources with or without Lyα. Since it is easier to observe Lyα for galaxies with higher SFR (at a given fesc), this already indicates that fesc is higher for galaxies with low SFR. Note that the SFR for the AGN is likely to be overestimated as AGN activity also contributes to the Hα flux.

Figure 6.

SFR(Hα) versus stellar mass for the observed HAEs. We obtained the SFR from dust-corrected Hα and stellar mass from SED fitting (see Section 3.5). We show the position of sources with and without Lyα, AGN with Lyα and AGN without Lyα. There is no obvious difference in the SFRs or stellar masses between sources with or without Lyα. Since it is easier to observe Lyα for galaxies with higher SFR (at a given fesc), this already indicates that fesc is higher for galaxies with low SFR. Note that the SFR for the AGN is likely to be overestimated as AGN activity also contributes to the Hα flux.

Fig. 5 shows the distribution of Hα fluxes of our observed sample, of the AGN and also of the star-forming sources directly detected in Lyα. Whether a source is detected in Lyα does not clearly correlate with Hα flux, such that even very faint HAEs are detected. These very faint sources generally have high fesc, although we note that it is possible that these sources are detected at a redshift where the transmission in the Hα filter is low. In the remainder of the paper, we will use the sample of 17 dual Hα-LAEs for direct measurements of fesc, and use upper limits for the other 471 HAEs.

### Lyα properties of dual Lyα-HAEs

After excluding the X-ray AGN, we find 12 robust dual Lyα-HAEs, which would translate in 2.5 per cent of our SFGs being detected in Lyα down to our flux limit. However, if we only select the subset of HAEs with the deepest Lyα observations (194 SFGs, of which eight have Lyα), we find a fraction of 4.1 per cent. This is comparable to Oteo et al. (2015), who found a fraction of 4.5 per cent and lower than the 10.9 per cent of Hayes et al. (2010), whose Lyα observations are a factor of 6 deeper and Hα sample consists of fainter sources (by a factor of 7), but covers a volume which is ∼80 times smaller than our survey.

Comparing our 12 Hα-LAEs to HAEs without Lyα, we find that there are no clear differences in the SFR-Mstar plane (Fig. 6). Since it is easier to observe Lyα for galaxies with higher SFR (at a given fesc), this already indicates that fesc is higher for galaxies with low SFR. We show the stellar masses, observed Hα and Lyα fluxes and dust attenuations of the Lyα detected sources in Table 2.

The median, dust-corrected Hα luminosity for the 12 SFGs with Lyα is 4.5 × 1042 erg s−1(corresponding to a SFR of 20 M yr−1), and median stellar mass of 1.8 × 1010 M, such that they have specific SFRs which are typical to the sample of HAEs (see Fig. 6). The median Lyα luminosity is 2.8 × 1042 erg s−1 and the Lyα escape fraction for these galaxies ranges from 1.9 ± 0.3 to 30 ± 3.8 per cent (although we note this does not take the uncertainty due to the filter transmission profiles into account, nor a statistical correction.). The EW0(Lyα) ranges from 10 to 244 Å. If we apply an EW0(Lyα) cut of >25 Å (similar to the selection of LAEs at high redshift; e.g. Matthee et al. 2014 and references therein), seven out of 194 star-forming HAEs with deepest Lyα observations are recovered as LAEs, with luminosities ∼2–8 × 1042 erg s−1.

#### Morphology

In Fig. 7, we compare the Lyα SB with rest-frame UV, R and Hα from HST ACS F814W (Koekemoer et al. 2007), K and NBK (continuum corrected with K), respectively. In order to be comparable, the PSF of the HST images is matched to that of the NB392 imaging on the INT, using a convolution with a Gaussian kernel. As the PSF of our INT imaging is 1.8–2.0 arcsec, this is a major limitation. However, for the most significantly detected sources (in Lyα), it is still possible to study differences qualitatively. Even though there is ground-based I-band data available, we use HST data because those are deeper.

Figure 7.

Lyα morphologies for the 10 HAEs with most significant Lyα detections. All images (including HST ACS F814W) were smoothed to the PSF of the NB392 image and are centred on the position of peak Hα emission. The top row shows AGN, and the bottom row SFGs. The sources are ordered by Hα flux, decreasing from left to right. Each panel shows a 100 × 100 kpc K band (which traces rest-frame R, and thus roughly stellar mass) thumbnail with contours of rest-frame Lyα (blue), UV (green, from ACS F814W) and Hα (red). The Lyα image was obtained by subtracting the continuum PSF matched U band from the NB392 image and correcting using the differential PSF image. The contours show 3, 4, 5 and 6 σ levels in each respective filter. For Lyα, the 3σ contour corresponds to a SB ranging from 1.8 to 17 (median 2.0) × 10−18 erg s−1 cm−2 arcsec−2. The UV data are typically 3 mag deeper than the Lyα data. Among the AGN, IDs 7801 and 1139 show evidence for extended Lyα emission. For the SFGs there is little convincing evidence for extended emission at the reached SB limit. The Lyα emission of source 1993 appears to be offset from the peak UV emission. The faintest HAEs have no 3σ contour because of smoothing with the PSF of the Lyα image.

Figure 7.

Lyα morphologies for the 10 HAEs with most significant Lyα detections. All images (including HST ACS F814W) were smoothed to the PSF of the NB392 image and are centred on the position of peak Hα emission. The top row shows AGN, and the bottom row SFGs. The sources are ordered by Hα flux, decreasing from left to right. Each panel shows a 100 × 100 kpc K band (which traces rest-frame R, and thus roughly stellar mass) thumbnail with contours of rest-frame Lyα (blue), UV (green, from ACS F814W) and Hα (red). The Lyα image was obtained by subtracting the continuum PSF matched U band from the NB392 image and correcting using the differential PSF image. The contours show 3, 4, 5 and 6 σ levels in each respective filter. For Lyα, the 3σ contour corresponds to a SB ranging from 1.8 to 17 (median 2.0) × 10−18 erg s−1 cm−2 arcsec−2. The UV data are typically 3 mag deeper than the Lyα data. Among the AGN, IDs 7801 and 1139 show evidence for extended Lyα emission. For the SFGs there is little convincing evidence for extended emission at the reached SB limit. The Lyα emission of source 1993 appears to be offset from the peak UV emission. The faintest HAEs have no 3σ contour because of smoothing with the PSF of the Lyα image.

The sources with IDs 2741, 78011073, 9274 and 1139 (see Table 2 for more information), shown in the first row, are all AGN with mostly symmetrical Lyα morphology. Compared to the UV, IDs 7801 and 1139 show evidence for extended Lyα emission, while this is more evident when Lyα is compared to Hα. Note that the Lyα image is typically the shallowest, and that the outer contours of Lyα therefore typically represent a higher fraction of the peak flux than the UV contours.

The sources in the second row are undetected in the X-ray, and therefore classed as SFGs. 1057 is relatively massive (Mstar = 1010.6 M), and has fesc of 10.8 ± 1.6 per cent. 7693 has an intermediate mass of ∼109.5 M and escape fraction of 5.5 ± 0.8 per cent. From comparison with the Hα image, it can be seen that Lyα preferentially escapes offset to the south from the galaxy centre, which might be indicative of an outflow.

The sources with IDs 2600, 1993 and 4861, in the last three columns of Fig. 7, are the faintest HAEs for which we detect Lyα directly, such that there are no 3σ Hα contours due to smoothing the image with the PSF of the Lyα image. These HAEs have (possibly) little dust, blue UV slopes and Lyα EW0 > 30 Å, such that they would be selected as LAE. The masses, SFRs and blue UV slopes are consistent with results from typical LAEs (e.g. Nilsson et al. 2009; Ono et al. 2010), and similar to simulated LAEs (e.g. Garel et al. 2012, 2015). The Lyα emission for ID 1993 and 4861 appears to be offset from the peak UV emission. This indicates that slit spectroscopy of UV or Hα selected galaxies might miss significant parts of Lyα. Note that we look at the stacked UV, Lyα, Hα and morphologies of these 12 SFGs and the full sample of SFGs in Section 7.

### [O ii] and [O iii] emission lines

In addition to NBK imaging, the HiZELS survey consists of NBJ and NBH imaging in the same fields. These filters are designed such that they also cover the redshifted [O ii] (in NBJ) and [O iii] (in NBH) emission at z = 2.23, similar to e.g. Nakajima et al. (2012) and Sobral et al. (2012). Out of the 488 HAEs, 23 and 70 galaxies are detected in [O ii] and [O iii], respectively. Two out of the nine AGN are detected as line-emitters in both NBs, two as [O iii] emitters, and two AGN as [O ii] emitters. The three remaining AGN all have a spectroscopic redshift at z = 2.21–2.23 (Sobral et al. 2016). This means that all our AGN are either spectroscopically confirmed or have highly accurate photometric redshifts, with emission lines in at least two NBs.

There are two star-forming HAEs with Lyα and [O iii] (ID 1993 and 2600; see Fig. 7 and Table 2). Compared to all HAEs with detection in [O iii] and similar limits on the Lyα flux and fesc, these galaxies have the lowest mass and SFR (∼108.6 M, SFR < 8 M yr−1), most extreme UV slopes (either bluest, ID 2600, or reddest, ID 1993). They have fesc ∼ 10 per cent and [O iii] EW0 ∼ 100 Å, resembling local Green pea galaxies (e.g. Henry et al. 2015). Their dust corrections are uncertain, because the two methods described in Section 3.4.3 give a factor of 3–5 difference depending on the method. Compared to the other Lyα detected HAEs, these galaxies have the highest Lyα EWs. ID 2600 has very high Hα EW0 of ∼960 Å, similar to the most extreme emission line galaxies (e.g. Amorín et al. 2015), while the Lyα EW0 is 80 Å. ID 1993, however, has a relatively low Hα EW0 of 70 Å, and Lyα EW0 of 68 Å. We note that if we use the Garn & Best (2010) dust correction, these two galaxies would have a factor of 2–3 higher fesc.

## Correlations between Lyα escape and galaxy properties

We use our sample of dual Lyα-HAEs and upper limits for the other 468 star-forming HAEs to search for potential correlations between galaxy properties and the Lyα escape fraction, shown in Figs 8 and 9. Our sample is compared with the Hα selected sample of Oteo et al. (2015) and the spectroscopically, Lyα-selected sample of Song et al. (2014). The main difference with our sample is that it has deeper Lyα observations than Oteo et al. (2015) and spans a wider range of galaxy properties. The major difference in respect to Song et al. (2014) is that their sample is selected on strong Lyα, and therefore biased towards sources with high fesc.

Figure 8.

Lyα escape fraction versus SFR and stellar mass for galaxies without AGN for individual sources (left panels) and stacks (right panels). The green circles show our directly detected Hα-LAEs; grey triangles highlight upper limits (green triangles have a UV detection in NB392) and our X-ray identified AGN with Lyα are shown in red diamonds. We also add the Hα selected sample from Oteo et al. (2015) in orange stars and the Lyα selected sample from Song et al. (2014) as blue pentagons. Our survey clearly extends the probed parameter space in galaxy properties. Stacked values are shown for two different radial distances to the centre (corresponding to 3 arcsec diameter apertures – diamonds, and 6 arcsec diameter apertures – pentagons). The typical measurement uncertainty in stellar mass is indicated in the bottom left panel. Top row: the left panel shows the SFR obtained from Hα versus fesc of individual sources. Although a correlation is expected by definition, it can be seen that on average, galaxies with a higher SFR have a lower escape fraction. The grey region in the right panel shows the SFRs typical for galaxies in the sample from Hayes et al. (2010), who inferred fesc = 5.3 ± 3.8 per cent. Bottom row: escape fraction versus stellar mass. While fesc can be relatively high for low-mass galaxies, the stacked results indicate that fesc might decrease weakly with increasing stellar mass. The large difference in fesc between some massive individual sources and the stacked values indicates that there is likely significant scatter in the values of fesc at this mass range.

Figure 8.

Lyα escape fraction versus SFR and stellar mass for galaxies without AGN for individual sources (left panels) and stacks (right panels). The green circles show our directly detected Hα-LAEs; grey triangles highlight upper limits (green triangles have a UV detection in NB392) and our X-ray identified AGN with Lyα are shown in red diamonds. We also add the Hα selected sample from Oteo et al. (2015) in orange stars and the Lyα selected sample from Song et al. (2014) as blue pentagons. Our survey clearly extends the probed parameter space in galaxy properties. Stacked values are shown for two different radial distances to the centre (corresponding to 3 arcsec diameter apertures – diamonds, and 6 arcsec diameter apertures – pentagons). The typical measurement uncertainty in stellar mass is indicated in the bottom left panel. Top row: the left panel shows the SFR obtained from Hα versus fesc of individual sources. Although a correlation is expected by definition, it can be seen that on average, galaxies with a higher SFR have a lower escape fraction. The grey region in the right panel shows the SFRs typical for galaxies in the sample from Hayes et al. (2010), who inferred fesc = 5.3 ± 3.8 per cent. Bottom row: escape fraction versus stellar mass. While fesc can be relatively high for low-mass galaxies, the stacked results indicate that fesc might decrease weakly with increasing stellar mass. The large difference in fesc between some massive individual sources and the stacked values indicates that there is likely significant scatter in the values of fesc at this mass range.

Figure 9.

Correlations between fesc and dust extinction, rest-frame BV colour and β, with symbols as defined in Fig. 8. We indicate the maximum typical uncertainties on extinction and β in the left panels. Top row: the left panel shows that fesc is anti-correlated with the dust extinction, although there is significant scatter. This means that at fixed aperture, Lyα preferentially escapes galaxies with low dust content. For comparison, we show the relation at z ∼ 2–3 from Hayes et al. (2011) and z ∼ 0.3 from Atek et al. (2014). Our best-fitted relation (to detections and upper limits; see text) resembles more than that of local galaxies. The grey band shows the 1σ error of the normalization of the fit. Right: stacked values of escape fraction versus dust extinction. A similar trend is seen as for individual sources. However, because stacks are not biased towards high fesc values, the normalization is lower. Bottom row: escape fraction versus β for individual sources (left), which confirms that there is a bimodal relation between fesc and galaxy UV colour, such that either very blue or very red galaxies emit significant Lyα. As seen in stacks in the right panel, this bimodal trend is most clear at the largest apertures.

Figure 9.

Correlations between fesc and dust extinction, rest-frame BV colour and β, with symbols as defined in Fig. 8. We indicate the maximum typical uncertainties on extinction and β in the left panels. Top row: the left panel shows that fesc is anti-correlated with the dust extinction, although there is significant scatter. This means that at fixed aperture, Lyα preferentially escapes galaxies with low dust content. For comparison, we show the relation at z ∼ 2–3 from Hayes et al. (2011) and z ∼ 0.3 from Atek et al. (2014). Our best-fitted relation (to detections and upper limits; see text) resembles more than that of local galaxies. The grey band shows the 1σ error of the normalization of the fit. Right: stacked values of escape fraction versus dust extinction. A similar trend is seen as for individual sources. However, because stacks are not biased towards high fesc values, the normalization is lower. Bottom row: escape fraction versus β for individual sources (left), which confirms that there is a bimodal relation between fesc and galaxy UV colour, such that either very blue or very red galaxies emit significant Lyα. As seen in stacks in the right panel, this bimodal trend is most clear at the largest apertures.

In addition to studying correlations for individual sources, we also stack different subsets of galaxies. As described in Section 4, we will limit ourselves to the HAEs with the deepest NB392 observations. We divide our sample by SFR(Hα), stellar mass, dust extinction and UV slope and ensure that our results do not depend on our particular choice of bin limits and width by perturbing both significantly. The benefit from studying correlations in bins of galaxies is that the results are less dependent on the systemic uncertainties in the dust corrections. While the systematic difference in dust corrections can be up to a factor of 5 for individual sources, the differences are much smaller over a statistical sample (compare for example Tables 1 and 2).

### Varying SFR(Hα) and mass

In the top row of Fig. 8 we show that fesc is anti-correlated with SFR. This is seen both in the individual sources and in stacks. As both fesc and SFR involve the Hα measurement, we naively expected to find an anti-correlation between the two. Yet, the quantitative behaviour of this trend is not a priori trivial due to complex Lyα radiative transport. The qualitative trend is not strongly sensitive to the dust correction method applied. By fitting a linear relation to our stacked values (in logarithmic space), we obtain

(5)
\begin{eqnarray} {\rm log}_{10}({f}_{\rm esc}/({per cent}))= 1.34^{+0.12}_{-0.12} - 1.32^{+0.10}_{-0.10}\, {\rm log}_{10}\nonumber\\ (\rm SFR/(M_{{\odot }}\rm yr^{-1})) \end{eqnarray}
for a 12 kpc aperture radius, with $$\chi ^2_{\rm red} = 1.96$$. At a radius of 24 kpc, the normalization is $$2.43^{+0.15}_{-0.15}$$ and slope $$-1.65^{+0.08}_{-0.09}$$, with a $$\chi ^2_{\rm red} = 1.27$$. This means that there is an ∼13σ anti-correlation between the escape fraction and the SFR (as a slope of zero is rejected at 13σ). According to this relation, a galaxy with a SFR of 4 M yr−1 has a typical Lyα escape fraction of ∼3.5 ± 1.8 per cent (in 12 kpc apertures). We note that this relation probably turns over at lower SFR than we probe (as fesc would otherwise reach values >100 per cent). The top-right panel of Fig. 8 also shows that the Lyα escape fraction is higher at larger radii at all SFRs, although the errors become a bit larger.

It can be seen that our directly detected sample occupies the region in parameter space which has highest escape fraction at fixed SFR. Compared to the Lyα selected sample from Song et al. (2014), we find a stronger anti-correlation between fesc and SFR at a lower normalization. This is because Lyα selected sources are biased towards high values of fesc.

The bottom row of panels in Fig. 8 shows how fesc is related to stellar mass. While the stacked values show a weak anti-correlation (although at very low significance), there is no evidence for a trend between fesc and stellar mass for the individually detected sources even though our sample clearly extends the probed dynamic range up to higher masses. As massive galaxies would naively be expected to have a lower fesc, since they tend to have a higher SFR and are dustier (e.g. Ibar et al. 2013), this means that the Lyα escape fraction is not only determined by the scale of a galaxy and it is likely that more subtle processes such as dust and gas dynamics play an important role. Interestingly, the individual detected sources with high stellar mass are at much higher fesc than the stacked values, indicating that there is significant scatter in fesc in this mass range.

### Dust extinction and UV slope

In the top row of Fig. 9 we show how fesc is related to dust extinction. For individual sources and for stacks, we find that there is an anti-correlation between fesc and dust extinction, although this relation contains significant scatter, caused in part by errors in our measurement of fesc and E(BV). With our data, we are able to extend the previous found relation from Hayes et al. (2011) to lower escape fractions and higher dust extinctions. We fit the following linear relation:

(6)
$${f}_{\rm esc} = C \times 10^{-0.4 {\rm E}(B-V) k_{\rm Ly\alpha }}.$$
This fit is performed by simulating a large grid of normalizations and slopes and computing the χ2 for each combination of normalization and slope in log(fesc) − E(BV) space. Upper limits are taken into account by assigning them an fesc of 0.01 per cent and using their upper limit as error. Since a fit to individual galaxies is mostly determined by the directly detected dual Hα-LAEs, which are biased towards high fesc, we also fit to the stacked values.

By minimizing the χ2 for individual galaxies, we find $$C = 0.17^{+0.15}_{-0.09}$$ and $$k_{\rm Ly\alpha } =5.60^{+3.45}_{-2.90}$$, such that a galaxy with E(BV) = 0 has an escape fraction of 17 per cent. This is lower in normalization (although the errors are significant), and significantly shallower in slope than the fit from Hayes et al. (2011), who find C is 0.445 and kLyα = 13.8. The normalization and slope are more similar to the z = 0.3 result from Atek et al. (2014) (C = 0.22, kLyα = 6.67). A possible explanation could be that Hayes et al. (2011) miss dusty galaxies, such that they infer a steeper slope, which would be consistent with the discussion in Oteo et al. (2015). We discuss this further in Section 8.2. For stacks, we find $$C = 0.03^{+0.01}_{-0.01}$$ and $$k_{\rm Ly\alpha } =10.71^{+0.89}_{-1.01}$$ for 12 kpc apertures and $$C = 0.08^{+0.02}_{-0.01}$$ and $$k_{\rm Ly\alpha } =7.64^{+1.38}_{-1.36}$$ for 24 kpc apertures. Our fit to stacked data is less biased towards high values of fesc and therefore at a lower normalization. The slope is slightly higher, although still not as high as the slope inferred by Hayes et al. (2011). Similar as seen for the stacks in bins of stellar mass, the individually detected galaxies at highest dust attenuations have much higher fesc than the median stack. This means that there is a lot of scatter in the values of fesc at the highest dust attenuations.

We furthermore note that part of the correlation between fesc and E(BV) is expected because there is a dust correction in the Hα flux, and thus in fesc. The fact that the correlation is rather weak (the slope is inconsistent with being zero at only 1.9σ and as there is significant scatter) indicates that dust is not the only regulator of Lyα escape, a result already found by Atek et al. (2008) at low redshift. We also note that the trend between fesc and E(BV) becomes somewhat bimodal when using the Garn & Best (2010) dust correction, meaning that there are galaxies with high dust attenuation and high fesc, which is virtually impossible with a Calzetti et al. (2000) dust correction. Galaxies with low E(BV) can in this case also have a lower fesc, leading to a flattening of the relation.

The bottom row of panels in Fig. 9 shows how fesc is related with the UV slope β. As the UV slope is also a tracer of dust attenuation, we also find a tentative anti-correlation with escape fraction for galaxies with β < 0, but for stacks and individual sources, particularly when the HAEs from Oteo et al. (2015) are included. Surprisingly, the trend between fesc and β seems to reverse for redder galaxies, leading to a bimodal relation which is particularly seen in measurements of fesc with 3 arcsec apertures. The (maximum) typical error of β is indicated in the bottom left panel of Fig. 9 and is sufficiently small to exclude measurement errors of β as the source of the observed bimodal behaviour. The bimodal behaviour is also seen when correcting for dust using the Garn & Best (2010) prescription.

## EXTENDED EMISSION

As Lyα is a resonant emission line, scattering due to neutral hydrogen leads to a diffusion process similar to a random walk, which results in a lower SB. Therefore, it is likely that in the presence of extended Hi, Lyα emission will be more extended than the UV. Although our relatively large PSF FWHM and limited depth in NB392 (∼1.8 arcsec and ∼24 AB magnitude, respectively) limit the study of extended Lyα emission at low SB, we can test how our measured escape fraction depends on the chosen aperture size. We do this by analysing the stacked images for the (biased) sample of 12 dual Hα-Lyα emitting SFGs, and for the stack of the 265 SFGs (see Section 5 and Table 1). We measure both Hα, Lyα and the rest-frame UV in apertures ranging from 2 to 10 arcsec in diameter.

As seen in Fig. 10 (and illustrated by Fig. 11), Lyα is significantly more extended than the UV (traced by convolved HST F814W imaging) and Hα for both stacks. The stack of all SFGs (right panel of Fig. 11) has extended Lyα emission up to ∼20 kpc distance from the centre, at 3σ. At this significance, the stack of the 12 dual Hα-LAEs (left panel of Fig. 11) is extended up to 30 kpc, and is clearly more extended than the aperture that we used for the sources individually.

Figure 10.

Fraction of the total retrieved flux (defined as the maximum flux within 5 arcsec radius) as a function of aperture for the UV, Hα and Lyα and stack of direct detected Hα-LAEs (left) and stack of all SFGs (right). It can be seen that while Hα and the UV (from PSF convolved HST F814W imaging) quickly reaches the total flux in both cases, the Lyα flux continues to increase. For the direct detected sources, the Lyα flux increases somewhat more rapidly because sources were selected on bright central Lyα. These sources also have more compact Hα than UV emission. For the stack of all SFGs, Lyα continues to increase (although the errors are significant), such that deeper observations are required to test whether we have captured the full extent of Lyα.

Figure 10.

Fraction of the total retrieved flux (defined as the maximum flux within 5 arcsec radius) as a function of aperture for the UV, Hα and Lyα and stack of direct detected Hα-LAEs (left) and stack of all SFGs (right). It can be seen that while Hα and the UV (from PSF convolved HST F814W imaging) quickly reaches the total flux in both cases, the Lyα flux continues to increase. For the direct detected sources, the Lyα flux increases somewhat more rapidly because sources were selected on bright central Lyα. These sources also have more compact Hα than UV emission. For the stack of all SFGs, Lyα continues to increase (although the errors are significant), such that deeper observations are required to test whether we have captured the full extent of Lyα.

Figure 11.

Thumbnails of the stacks of the 12 Hα-LAEs (left) and full sample of SFGs (right). The background images show the K band (where we removed the contribution from Hα using NBK), which corresponds to rest-frame R. The Lyα emission is shown in blue contours, Hα emission is shown in red and UV (from observed F814W, corresponding to ∼2000 Å rest frame) in green. In the left panel we also indicate the 3 arcsec diameter aperture used for measurements of individual sources in black dashed lines. The PSF FWHM is indicated with a white dashed circle in both panels. The contour levels are normalized to the peak flux in each band. The outer Lyα contour represents 7.5σ for the left panel, and 3σ for the right panel. This corresponds to 0.37 and 0.52 of the peak flux, respectively. The 1σ SB limit in the two panels is 9.0 and 2.5 × 10−19 erg s−1 cm−2 arcsec−2. Other contours correspond to a fraction of 0.5 and 0.75 of the peak flux in the left panel, and 0.6, 0.75 and 0.9 in the right panel. In both panels, it can be seen that Hα traces the UV light very well. Lyα is more extended than Hα and the UV for both the (biased) stack of direct detect Hα-Lyα and the full stack of SFGs, indicative of scattering. Lyα extends up to ∼20 kpc distance from the centre at the corresponding significances. The stack of the 12 Hα-LAEs is extended up to 30 kpc at 3σ significance.

Figure 11.

Thumbnails of the stacks of the 12 Hα-LAEs (left) and full sample of SFGs (right). The background images show the K band (where we removed the contribution from Hα using NBK), which corresponds to rest-frame R. The Lyα emission is shown in blue contours, Hα emission is shown in red and UV (from observed F814W, corresponding to ∼2000 Å rest frame) in green. In the left panel we also indicate the 3 arcsec diameter aperture used for measurements of individual sources in black dashed lines. The PSF FWHM is indicated with a white dashed circle in both panels. The contour levels are normalized to the peak flux in each band. The outer Lyα contour represents 7.5σ for the left panel, and 3σ for the right panel. This corresponds to 0.37 and 0.52 of the peak flux, respectively. The 1σ SB limit in the two panels is 9.0 and 2.5 × 10−19 erg s−1 cm−2 arcsec−2. Other contours correspond to a fraction of 0.5 and 0.75 of the peak flux in the left panel, and 0.6, 0.75 and 0.9 in the right panel. In both panels, it can be seen that Hα traces the UV light very well. Lyα is more extended than Hα and the UV for both the (biased) stack of direct detect Hα-Lyα and the full stack of SFGs, indicative of scattering. Lyα extends up to ∼20 kpc distance from the centre at the corresponding significances. The stack of the 12 Hα-LAEs is extended up to 30 kpc at 3σ significance.

The growth curves for Hα and the UV are similar, quickly growing to the maximum at ∼20 kpc. Lyα, however, continues to increase. The Lyα flux for the stack of the sample of 12 Hα-LAEs peaks at ∼4 arcsec radial apertures, as further increase in the fraction of recovered flux with increasing aperture is within the errors. The fesc within this radius, ∼30 kpc, is 14.2 ± 1.9 per cent. At a radius of 1.5 arcsec (similar to the aperture used for individual sources), we measure a stacked fesc of 7.7 ± 0.9 per cent. At this aperture, only ∼50 per cent of the maximum observed Lyα flux is retrieved. The Lyα flux for the stack of all SFGs also continues to increase up to at least 30 kpc. By fitting a linear relation between the fraction of the total recovered flux and radius, we find that r90 (the radius at which 90 per cent of the flux is retrieved) for our stack of directly detected Hα-Lyα is 31.3 ± 1.5 kpc. For the stack of all SFGs, r90 = 36.0 ± 3.8 kpc. For Hα, we find values of r90 = 19.3 ± 1.6 kpc and r90 = 21.0 ± 0.5 kpc, respectively.

We show the SB profile of the full stack of SFGs in the left panel of Fig. 12, where we scaled Hα such that it has a similar SB as Lyα at 5 arcsec radius, corresponding to an escape fraction of 2 per cent (see also the right panel). The Hα profile follows an exponential (as the y-axis is logarithmic), decreasing with increasing radius. While the Lyα profile seems to be more complex, we note that the errors due to the different PSF shapes of broad-band and NB are important, particularly in the centre part. Towards higher radii (>30 kpc), the Lyα signal from the stack is significantly above the errors due to differences in the PSF, and is thus real. We find no evidence that the integrated Lyα flux does not continue to grow up to 40 kpc distance, indicating that it can be extended up to larger radii if deeper SB limits are reached. This also means that we cannot yet directly infer the total fesc.

Figure 12.

Left: SB of our full stack of SFGs for Hα and Lyα and we also show the SB for the stack of the PSF control sample (Section 4.1), for which the stack of SFGs is corrected. The SB that we observe is attributed to different shapes of the PSF of the NB and broad-band. This SB difference is added to the error of the Lyα SB. We scaled the Hα flux such that it has the same SB as Lyα at a radius of 5 arcsec, corresponding to 41 kpc. Note that this scaling corresponds to an escape fraction of 2 per cent, discussed in Section 7.1. Right: escape fraction of the stack of all sources as a function of radial distance from the centre. We compare stacking sources with and without removal of AGN, and we also show the escape fraction if we do not correct Hα for dust. The escape fraction at large radii is slightly higher without the removal of AGN, because the AGN are typically bright in Lyα and the haloes around them have a larger scale.

Figure 12.

Left: SB of our full stack of SFGs for Hα and Lyα and we also show the SB for the stack of the PSF control sample (Section 4.1), for which the stack of SFGs is corrected. The SB that we observe is attributed to different shapes of the PSF of the NB and broad-band. This SB difference is added to the error of the Lyα SB. We scaled the Hα flux such that it has the same SB as Lyα at a radius of 5 arcsec, corresponding to 41 kpc. Note that this scaling corresponds to an escape fraction of 2 per cent, discussed in Section 7.1. Right: escape fraction of the stack of all sources as a function of radial distance from the centre. We compare stacking sources with and without removal of AGN, and we also show the escape fraction if we do not correct Hα for dust. The escape fraction at large radii is slightly higher without the removal of AGN, because the AGN are typically bright in Lyα and the haloes around them have a larger scale.

Comparing the Hα and Lyα profiles results in an increasing Lyα escape fraction with increasing aperture (see the right panel of Fig. 12). The fesc increases from 0.3 ± 0.05 per cent at 12 kpc distance to 1.6 ± 0.5 per cent at ∼30 kpc. Note that without dust correction fesc is roughly a factor of 2 higher. At the radius of 30 kpc, the Lyα SB is ∼6 × 10−19 erg s−1 cm−2 arcsec−2 (see also the left panel), such that the extended emission is detected at ∼2.4σ (at 2σ confidence level, extended emission is seen up to ∼40 kpc). As can be seen in the right panel of Fig. 10 and illustrated in the right panel of Fig. 11, Lyα is extended up to ∼20 kpc at 3σ confidence level. At 2σ significance, it is extended up to ∼30 kpc. This means that aperture-based measurements (including slit spectroscopy) might miss parts of Lyα emission, and that IFU's or specially designed NB filters are more suited.

## DISCUSSION

### A consensus on the value of fesc

To date, a number of papers have been published on measuring the Lyα escape fraction with different selections of galaxies and methods. Typically, Lyα selected galaxies at z > 2 have resulted in high escape fractions of ∼30 per cent (e.g. Blanc et al. 2011; Nakajima et al. 2012; Wardlow et al. 2014; Kusakabe et al. 2015; Trainor et al. 2015), even though techniques to estimate the intrinsic Lyα production range from UV to dual NB to FIR and spectroscopy. However, a Lyα selected sample of galaxies is not representative for the full (star-forming) galaxy population and estimates of fesc are biased towards high values (as e.g. seen in Fig. 8).

With UV or rest-frame optical emission line selected galaxies, the typical fesc at z = 2–3 is around fesc ∼ 3–5 per cent (e.g. Hayes et al. 2010; Kornei et al. 2010; Ciardullo et al. 2014), but fesc is found to increase with increasing redshift, up to ∼30 per cent at z = 5.7 (e.g. Hayes et al. 2011). Our measured median value of 1.6 ± 0.5 per cent of HAEs is lower than the value in these papers in the literature. However, we note that our measurement is the first which is independent of assumptions on the shape of the luminosity function and integration limits. Furthermore, as we will show now, the results are fully consistent with literature results when we account for different selections of galaxies and the different parameter spaces probed by different surveys.

Most samples have either been UV selected or selected with emission lines bluer than Hα (e.g. [O iii]; Ciardullo et al. 2014), such that they might miss a population of dusty galaxies. This might cause a bias towards high fesc values, as we show here that dustier galaxies typically have a lower fesc (see Fig. 9, but note that we also observe significant fesc for some dusty galaxies). On the contrary, Cassata et al. (2015) report a low fraction of UV selected galaxies with Lyα emission at 2 < z < 6, corresponding to fesc < 1 per cent, concluding that the bulk of the Lyα luminosity density is coming from fainter galaxies than the typically UV bright galaxies that were targeted. This is consistent with our results, as these galaxies have a typical SFR of ∼50–100 M yr−1 (Tasca et al. 2015), and we show that the typical escape fraction at these SFRs is very low (<0.5 per cent).

The major difference between our survey and the matched Hα-Lyα survey at z = 2.2 from Hayes et al. (2010) is that our HAEs probe a larger range in SFRs and stellar masses, because our observations probe a much larger volume (∼80 times larger, yet with shallower depth). We have shown in the top row of panel of Fig. 8 that the SFR anti-correlates with the escape fraction. If we compute SFRs from the HAEs from Hayes et al. (2010) as described in Section 3.5, we find that their SFRs are between 0.4 and 4.4 M yr−1, illustrated as the grey region in the top-right panel of Fig. 8. The difference between the inferred escape fraction of 5.3 ± 3.8 per cent by Hayes et al. (2010) and our stack of SFGs can be fully explained by the different parts of parameter space covered, and the anti-correlation between fesc and SFR. This also means that the bulk of Lyα luminosity density at z = 2.23 is not emitted by SFGs of >3 M yr−1 (roughly the lower limit of the Hα survey), but by much fainter galaxies. We note that Konno et al. (2015) infer a global fesc of ∼1.5 per cent based on more recent measurements of UV and Lyα luminosity functions.

### Dependence on galaxy properties

In Section 6 we have shown that fesc generally increases with decreasing SFR and dust attenuation (albeit with significant scatter), and that it shows bimodal behaviour with UV slope. The trends between fesc and SFR and dust attenuation are also seen in the local Universe (Hayes et al. 2014). In this subsection we discuss how our inferred trends compare to previous results and speculate on the origin of bimodal relations.

#### Relation with dust attenuation

Other surveys have also established a correlation between fesc and dust attenuation, both in the local Universe (Atek et al. 2008, 2014) and at high redshift (Verhamme et al. 2008; Hayes et al. 2010; Kornei et al. 2010; Blanc et al. 2011; Cassata et al. 2015). This trend is in line with the significantly lower observed fesc in the low-redshift Universe (e.g. ∼1 per cent at z = 0.3–1; Deharveng et al. 2008; Cowie, Barger & Hu 2010; Wold, Barger & Cowie 2014), as the dust content of galaxies increases with cosmic time (cf. Hayes et al. 2011).

We find a weak trend between fesc and dust attenuation, both for stacks and individual galaxies. As shown in the top row of panels of Fig. 9, our trends between fesc and E(BV) are significantly shallower than the trend inferred by Hayes et al. (2011), and resembles the trend seen in the local Universe (Atek et al. 2014). However, note that comparison is limited by systematic uncertainties and differences between the methods used to estimate fesc (e.g. Atek et al. 2014 use Balmer decrements to estimate the dust uncertainty, while our survey is dependent on SED fitting). While Hayes et al. (2010) only included Hα and Lyα selected galaxies in their analysis, Hayes et al. (2011) added the UV selected galaxies at z ∼ 3 from Kornei et al. (2010) in order to have a larger sample. However, the UV selected sample requires additional uncertainties on the intrinsic production of Lyα and is biased towards lower dust attenuations (even though their mass and SFR distribution are similar to our Hα selected sample).

We speculate that the prime origin of the discrepancy between our fitted trend and the trends from Hayes et al. (2010, 2011) is the combination that we probe more luminous galaxies with higher masses and SFRs (compared to their Hα selected sub-sample) and more dusty galaxies (compared to their UV selected sub-sample), and therefore find a lower normalization. As our parameter space also includes more dusty galaxies with significant fesc, the relation is flattened. We also note that we find that at a fixed dust attenuation, galaxies with a relatively higher escape fraction have a higher SFR, and are thus also more likely to not be present in the Hayes et al. (2011) sample (as galaxies with higher SFR are rarer), and leading to a flattening of the relation. At the highest dust attenuations, we observe individual galaxies with fesc over two orders of magnitude higher than the typical value for galaxies with similar dust attenuation, indicating that there is a lot of dispersion in the values of fesc.

It is interesting to note that we find that the E(BV) values of HAEs are not dependent on Lyα EW. The relation we would infer between fesc and E(BV) for HAEs with Lyα EW0 > 25 Å would be flatter than the relation inferred for our full sample of HAEs. This is because dusty HAEs with high fesc tend to have high Lyα EW0, which may be a sign of outflows. Note however that we do not include Lyα selected sources with faint Hα (and thus high fesc) in this analysis.

There still exist large systematic uncertainties on the method to correct for dust attenuation at high redshift. For example, if the nebular extinction is indeed stronger than the stellar extinction at high SFRs (e.g. Reddy et al. 2015), this may influence our observed trends in Figs 8 and 9. We use the Garn & Best (2010) prescription for dust attenuation to evaluate how such a differential dust correction affects our results, as this prescription is based on stellar mass and thus qualitatively similar to a dust attenuation which varies with SFR. On a source-by-source basis, the fesc can vary significantly (see e.g. the difference in attenuations in Table 2); however, the results for most stacks are very similar. The largest changes are seen for the stacks with either the highest or lowest SFRs, stellar masses and E(BV) values. We find that the relation between fesc and E(BV)star flattens somewhat, which is driven by a higher attenuation (and thus lower fesc) for sources in the bin with lowest E(BV) values. The relations between fesc and SFR and stellar mass would be steeper. Regardless of the method used to correct for dust, the trend we find between fesc and E(BV) is relatively weak and there is a large dispersion in the values of fesc. This indicates that dust attenuation alone is not the most important regulator of fesc.

#### Absence of relations and bimodality?

Additional evidence that dust is not the most important (or at least not the only) regulator of fesc is that there exists a significant population of red, dusty LAEs (e.g. Blanc et al. 2011; Guaita et al. 2011; Nilsson et al. 2011; Oteo et al. 2015; Taniguchi et al. 2015 and this survey), particularly in the presence of strong outflows (e.g. Atek et al. 2008; Martin et al. 2015). We find tentative evidence that there are bimodal trends between fesc and UV slope. This is supported by the observation that the two reddest galaxies in the sample (IDs 1057 and 1993) have a high fesc. These sources are surprisingly different: while 1993 has low mass and high dust content, 1057 is massive with little dust. These two sources are not atypical; in fact, almost all Lyα detected HAEs are either bluer or redder than the average HAE without Lyα detection (which have β ∼ −1). Apart from having high Lyα EW, we find no galaxy property which is related to having a red UV slope and high escape fraction. We note that this bimodal trend can only be seen for samples which include red, relatively massive objects. The trend is therefore still consistent with the results from Hagen et al. (2016), who find no statistical difference in the UV slopes of Lyα or [O iii] selected galaxies.

Fig. 8 shows that there is only a weak trend between fesc and mass for stacks, such that high mass galaxies have lower fesc. However, there are relatively massive galaxies which are detected with high fesc. This means that there is a lot of dispersion in the values of fesc at high stellar masses. This is also seen in the galaxies with highest dust attenuation. It is interesting to note that the mass scale at which we find higher fesc for individual sources coincides with the point where our Hα selection starts picking out galaxies below the main relation between stellar mass and SFR (see also Fig. 4). Because of this, and with observations of local galaxies with strong outflows and Lyα emission in mind, we speculate that this enhanced fesc might be related due to outflows of (dusty) gas. Outflows in turn redshift Lyα out of resonance wavelengths, which facilitate the escape of Lyα photons, as for example observed in local galaxies (e.g. Bik et al. 2015; Duval et al. 2015; Herenz et al. 2015; Rivera-Thorsen et al. 2015). This idea can be tested with spectroscopic follow-up.

### Extended Lyα emission

There is increasing evidence that SFGs are surrounded by a low SB Lyα halo, as indicated from stacks of UV selected galaxies (Steidel et al. 2011) and LAEs (Matsuda et al. 2012; Feldmeier et al. 2013; Momose et al. 2014). Lyα haloes have also been detected around faint individual Lyα selected galaxies, both locally and at high redshift (Hayes et al. 2013; Caminha et al. 2015; Patrício et al. 2016; Wisotzki et al. 2015). Moreover, luminous LAEs, without clear signs of AGN activity, can also be extended up to at least ∼16 kpc (e.g. Hayashino et al. 2004; Ouchi et al. 2013; Sobral et al. 2015b).

As Lyα emission is observed to be more extended than the UV, it is usually attributed to resonant scattering (e.g. Zheng et al. 2010; Dijkstra & Kramer 2012; Hayes et al. 2013), although simulations from Lake et al. (2015) indicate that collisional ionization due to cooling flows might also contribute. Very extended Lyα emission on scales of ∼100 kpc has been observed in Lyα blobs (e.g. Matsuda et al. 2004) and Lyα haloes around radio galaxies (e.g. Rottgering et al. 1995; Venemans et al. 2007; Swinbank et al. 2015), believed to be powered by either star formation or AGN activity (e.g. Geach et al. 2009; Overzier et al. 2013; Ao et al. 2015; Umehata et al. 2015).

We find that the average Hα selected SFG at z = 2.23 shows Lyα emission which is more extended than their Hα or the UV emission (see Figs 10–12). These galaxies are typically dustier than the sources from Steidel et al. (2011) (with our sample of HAEs having a median ALyα = 3.0, compared to ALyα = 1.12 from Steidel et al. 2011). The SB of the stack of all sources is therefore fainter (roughly a factor of 2, at a radial distance of 20 kpc) than the stack by Steidel et al. (2011), but we note that the stack of galaxies with E(BV) ∼ 0.1 (ALyα = 1.2) has similar Lyα SB values as the stack from Steidel et al. (2011). We also note that our stack of SFGs is fainter in the UV by 0.5 mag, but has a dust corrected SFR of ∼42 M yr−1, compared to the typical SFR of ∼34 M yr−1 of the sample from Steidel et al. (2011), because of the higher dust attenuation of our sample.

With the current depth, we cannot measure the maximum extent of Lyα for the stack of all SFGs or the stack of directly detected Hα-LAEs (see Fig. 10). By comparing the ratio of r90 of Lyα and Hα for both stacks, we find that Lyα radii are larger than Hα radii with a factor of 1.6–1.7. These values are comparable to those found for local Lyα analogues from Hayes et al. (2013), although that analysis used r20, which is undesirable for our sample due to the limitations from the PSF. For these local analogues, the average ratio between the Lyα and Hα radii is, however, slightly higher than the ratio for our stack of SFGs, which may indicate that we have not yet observed the full extent of Lyα.

For the stack of directly detected Hα-LAEs, we observed an escape fraction of fesc = 14.2 ± 1.9 per cent at a radius of 30 kpc, twice the fesc as the typical aperture used for photometry of Lyα selected galaxies. Assuming the Calzetti et al. (2000) dust law, we can estimate that attenuation of Lyα photons due to dust is a factor of 9 ± 3 (ALyα = 2.4 ± 0.3). Hence, this predicts that fesc before dust attenuation is 100/(9 ± 3) per cent. This indicates that dust attenuation is sufficient to explain that fesc is ∼15 per cent. This may well be the case for the majority of Lyα selected galaxies. However, we note that the sample of directly detected Hα-LAEs is obviously biased towards centrally peaked, high central fesc Lyα emission, as otherwise they would not have been directly detected in Lyα. It is therefore likely that the typical SFG has a larger Lyα halo with a less peaked central SB. The total inferred fesc is roughly a factor of 2 higher than fesc measured with 3 arcsec diameter apertures. This means that, if we naively assume that this factor of 2 is constant for all 12 SFGs used in the stack, the total escape fractions of our directly detected SFGs (see Table 2) range from 4 to 60 per cent.

## CONCLUSIONS

We have undertaken a panoramic matched Hα-Lyα NB survey to study the Lyα emission from a sample of 488 Hα selected SFGs at z = 2.23 and its dependence on radial scale and galaxy properties. Our conclusions are as follows.

• Out of the 488 observed HAEs, we detect 43 sources in the NB392 imaging. 17 of these have strong Lyα emission (of which five are X-ray AGN), and we measure the UV continuum for the other 26. We put limits on the escape fractions for these 26 and the other 445 undetected HAEs.

• The observed Lyα escape fraction for individual detected sources ranges from 2 to 30 per cent, and we find that these HAEs probe a wide range of star-forming systems: with masses from 3 × 108 M to 1011 M and dust attenuations E(BV) = 0 − 0.5. We particularly note that some massive, dusty galaxies are clearly visible in Lyα, while they have no evidence for AGN activity from X-ray observations.

• With matched NB observations in the J and H band from Sobral et al. (2013), we are also able to detect [O ii] and [O iii] emission lines for 23 and 70 galaxies, respectively. Two faint Lyα emitting galaxies are detected in [O iii]. Remarkably, these have relatively high escape fraction and have the lowest stellar mass and highest Lyα EW from our sample.

• While Lyα morphologies of individual X-ray selected AGN are typically circularly symmetric, tracing the rest-frame UV light, the Lyα morphologies of SFGs are more irregular. Lyα can appear to be offset from the rest-frame UV and Hα emission, and also more extended, although we detect this only significantly for AGN.

• Both for the individually detected sources and by stacking, we confirm existing trends between the Lyα escape fraction and SFR and dust attenuation and explore how they are extended to a larger range of parameter space – particularly more massive and dusty galaxies. Lyα escape increases for galaxies with lower SFR and low dust attenuation, albeit with significant scatter. The trend between fesc and dust attenuation resembles that observed in local galaxies. The escape fraction shows bimodal behaviour with UV slope: fesc increases for the bluest and reddest galaxies and is at a minimum UV slope β ∼ −0.5. At the highest masses and dust attenuations, we detect individual galaxies with fesc much higher than the typical values from stacking, indicating significant scatter in the values of fesc at the highest masses and dust attenuations.

• We interpret apparently contradictory reported numbers of the escape fraction by various studies using our observed trends with galaxy properties. Deeper surveys over a smaller area infer a higher escape fraction as the SFRs of their galaxies are lower in general. As these studies miss rarer galaxies with higher SFR and a higher dust attenuation, the inferred correlation between fesc and dust content steepens. Studies based on UV selected galaxies report higher numbers of fesc because the galaxies in their sample are biased towards galaxies with little dust, which tend to have higher fesc, except if their typical SFR is higher.

• The stack of our directly detected Hα-LAEs reveals that Lyα is more extended than the UV and Hα, with an extent of at least ∼30 kpc. For Lyα, the radius at which 90 per cent of the flux is retrieved is r90 = 31.3 ± 1.5 kpc, roughly 1.6 times larger than Hα. The median escape fraction for this biased sample increases up to 14.2 ± 1.9 per cent at least. We suggest that the missing Lyα photons can be fully explained by dust absorption. At the typical apertures used for slit spectroscopy (e.g. 1 arcsec), only half of the Lyα flux is recovered, even for the sample which is biased towards centrally peaked Lyα emission. This means that slit spectroscopy might miss a significant fraction of Lyα.

• By stacking our full sample of z = 2.23 SFGs, we find that Lyα extends up to at least 3 arcsec (with r90 = 36.0 ± 3.8 kpc) in radial distance at 3σ, corresponding to ∼20 kpc (30 kpc, 2σ), roughly twice the scale of star formation as traced by Hα and by the UV, with a Lyα SB of ∼6 × 10−19 erg s−1 cm−2 arcsec−2. As a result, the Lyα escape fraction of our full sample of HAEs continues to increase with radius, up to at least 1.6 ± 0.5 per cent at 30 kpc distance from the centre. With the current depth, it is not yet possible to constrain the maximum extent of Lyα.

Two important questions regarding the escape of Lyα remain. First, what is the total fesc for the typical SFG at z = 2.23? Deeper Lyα observations are required to fully answer this by measuring the full extent of Lyα. Secondly, properties which are not studied in this survey are likely equally or more important (e.g. the Hi covering fraction and geometry, e.g. Scarlata et al. 2009; Hi column density, Henry et al. 2015; outflows, Atek et al. 2008; ionization state, Hayes et al. 2014). From the significant scatter in the relations and from the bimodal behaviour in, for example, the relation between fesc and dust content and UV slope, it is clear that multiple mechanisms are likely together responsible for determining the Lyα escape fraction and its relation to extended Lyα emission.

It is crucial to improve the measurements of dust attenuation and to obtain precise measurements of redshifts and velocity offsets between Lyα and Hα, as these are responsible for the largest systematic errors which can go up to an order of magnitude when combined in unfortunate ways. However, the major downside is that this requires large, spectroscopic samples over a wide wavelength regime. Additional deeper Lyα NB observations are required in order to eliminate the biases caused by binning and stacking, and to be able to study the relation between galaxy properties and spatial extent of Lyα on a source-by-source basis.

We thank the anonymous referee for constructive comments and suggestions which have improved the quality of this work. JM acknowledges the support of a Huygens PhD fellowship from Leiden University. DS and JM acknowledge financial support from the Netherlands Organization for Scientific research (NWO) through a Veni fellowship, and DS from FCT through a FCT Investigator Starting Grant and Start-up Grant (IF/01154/2012/CP0189/CT0010) and from FCT grant PEst-OE/FIS/UI2751/2014. IO acknowledges support from the European Research Council (ERC) in the form of Advanced Investigator Programme, COSMICISM, 321302. HR acknowledges support from the ERC Advanced Investigator programme NewClusters 321271. IRS acknowledges support from STFC (ST/L00075X/1), the ERC Advanced Investigator programme DUSTYGAL 321334 and a Royal Society/Wolfson Merit Award. APA acknowledges support from the Fundaçao para a Ciência e para a Tecnologia (FCT) through the Fellowship SFRH/BD/52706/2014.

Based on observations made with the Isaac Newton Telescope (proposals 2013AN002, 2013BN008, 2014AC88, 2014AN002, 2014BN006, 2014BC118) operated on the island of La Palma by the Isaac Newton Group in the Spanish Observatorio del Roque de los Muchachos of the Instituto de Astrofísica de Canarias. We acknowledge the tremendous work that has been done by both COSMOS and UKIDSS UDS/SXDF teams in assembling such large, state-of-the-art multi-wavelength data sets over such wide areas, as those have been crucial for the results presented in this paper. The sample of HAEs is publicly available from Sobral et al. (2013).

We have benefited greatly from the publically available programming language Python, including the numpy, matplotlib, pyfits, scipy (Jones et al. 2001; Hunter 2007; Van Der Walt, Colbert & Varoquaux 2011) and astropy (Astropy Collaboration et al. 2013) packages, the imaging tools SExtractor, Swarp and Scamp (Bertin & Arnouts 1996; Bertin 2006, 2010) and the Topcat analysis program (Taylor 2005).

1
Note that we investigate the effect of different filter transmissions between Lyα and Hα as a function of redshift and the effect of systematic velocity offsets between the lines in Section 3.4.2.
2
Note that these sources are unlikely higher-redshift interlopers, as the NB392 filter is below the Lyman break at z > 3, such that a NB392 detection rules out that the source is a z ∼ 3.3 Hβ/[O iii] or z = 4.7 [O ii] emitter (e.g. Khostovan et al. 2015). Moreover, 13 of these have detected [O iii], [O ii] or both lines.

## REFERENCES

Amorín
R.
et al
2015
A&A

578
A105
Ao
Y.
et al
2015
A&A

581
A132
Astropy Collaboration
et al
2013
A&A

558
A33
Atek
H.
Kunth
D.
Hayes
M.
Östlin
G.
Mas-Hesse
J. M.
2008
A&A

488
491
Atek
H.
Kunth
D.
Schaerer
D.
Mas-Hesse
J. M.
Hayes
M.
Östlin
G.
Kneib
J.-P.
2014
A&A

561
A89
Bertin
E.
2006
Gabriel
C.
Arviset
C.
Ponz
D.
Enrique
S.
Astronomic Society of the Pacific Conference Series Vol. 351, Astronomical Data Analysis Software and Systems XV

Astron. Soc. Pac.
San Francisco
112
Bertin
E.
2010
SWarp: Resampling and Co-adding FITS Images Together, Astrophysics Source Code Library

(ascl:1010.068)
Bertin
E.
Arnouts
S.
1996
A&AS

117
393
Best
P.
et al
2013
Astrophys. Space Sci. Proc.

37
235
Bik
A.
Östlin
G.
Hayes
M.
A.
Melinder
J.
Amram
P.
2015
A&A

576
L13
Blanc
G. A.
et al
2011
ApJ

736
31
Boquien
M.
Buat
V.
Perret
V.
2014
A&A

571
A72
Bruzual
G.
Charlot
S.
2003
MNRAS

344
1000
Bunker
A. J.
Warren
S. J.
Hewett
P. C.
Clements
D. L.
1995
MNRAS

273
513
Calzetti
D.
Armus
L.
Bohlin
R. C.
Kinney
A. L.
Koornneef
J.
Storchi-Bergmann
T.
2000
ApJ

533
682
Caminha
G. B.
et al
2015
preprint (arXiv:1512.05655)
Cantalupo
S.
Arrigoni-Battaia
F.
J. X.
Hennawi
J. F.
P.
2014
Nature

506
63
Capak
P.
et al
2007
ApJS

172
99
Cassata
P.
et al
2015
A&A

573
A24
Chabrier
G.
2003
PASP

115
763
Ciardullo
R.
et al
2014
ApJ

796
64
Civano
F.
et al
2012
ApJS

201
30
Cowie
L. L.
Barger
A. J.
Hu
E. M.
2010
ApJ

711
928
Deharveng
J.-M.
et al
2008
ApJ

680
1072
Dijkstra
M.
2014
PASA

31
40
Dijkstra
M.
Kramer
R.
2012
MNRAS

424
1672
Dunlop
J. S.
McLure
R. J.
Robertson
B. E.
Ellis
R. S.
Stark
D. P.
Cirasuolo
M.
de Ravel
L.
2012
MNRAS

420
901
Duval
F.
et al
2015
preprint (arXiv:1512.00860)
Elvis
M.
et al
2009
ApJS

184
158
Erb
D. K.
et al
2014
ApJ

795
33
Feldmeier
J. J.
et al
2013
ApJ

776
75
Furusawa
H.
et al
2008
ApJS

176
1
Garel
T.
Blaizot
J.
Guiderdoni
B.
Schaerer
D.
Verhamme
A.
Hayes
M.
2012
MNRAS

422
310
Garel
T.
Blaizot
J.
Guiderdoni
B.
Michel-Dansac
L.
Hayes
M.
Verhamme
A.
2015
MNRAS

450
1279
Garn
T.
Best
P. N.
2010
MNRAS

409
421
Geach
J. E.
Smail
I.
Best
P. N.
Kurk
J.
Casali
M.
Ivison
R. J.
Coppin
K.
2008
MNRAS

388
1473
Geach
J. E.
et al
2009
ApJ

700
1
Gronwall
C.
et al
2007
ApJ

667
79
Guaita
L.
et al
2011
ApJ

733
114
Hagen
A.
et al
2016
ApJ

817
79
Hashimoto
T.
Ouchi
M.
Shimasaku
K.
Ono
Y.
Nakajima
K.
Rauch
M.
Lee
J.
Okamura
S.
2013
ApJ

765
70
Hayashino
T.
et al
2004
AJ

128
2073
Hayes
M.
2015
PASA

32
27
Hayes
M.
et al
2010
Nature

464
562
Hayes
M.
Schaerer
D.
Östlin
G.
Mas-Hesse
J. M.
Atek
H.
Kunth
D.
2011
ApJ

730
8
Hayes
M.
et al
2013
ApJ

765
L27
Hayes
M.
et al
2014
ApJ

782
6
Henry
A.
Scarlata
C.
Martin
C. L.
Erb
D.
2015
ApJ

809
19
Herenz
E. C.
et al
2015
preprint (arXiv:1511.05406)
Hummer
D. G.
Storey
P. J.
1987
MNRAS

224
801
Hunter
J. D.
2007
Comput. Sci. Eng.

9
90
Ibar
E.
et al
2013
MNRAS

434
3218
Ilbert
O.
et al
2009
ApJ

690
1236
Jones
E.
et al
2001
SciPy: Open source scientific tools for Python

Kashikawa
N.
et al
2012
ApJ

761
85
Khostovan
A. A.
Sobral
D.
Mobasher
B.
Best
P. N.
Smail
I.
Stott
J. P.
Hemmati
S.
Nayyeri
H.
2015
MNRAS

452
3948
Koekemoer
A. M.
et al
2007
ApJS

172
196
Konno
A.
Ouchi
M.
Nakajima
K.
Duval
F.
Kusakabe
H.
Ono
Y.
Shimasaku
K.
2015
preprint, (arXiv:1512.01854)
Kornei
K. A.
Shapley
A. E.
Erb
D. K.
Steidel
C. C.
Reddy
N. A.
Pettini
M.
Bogosavljević
M.
2010
ApJ

711
693
Kulas
K. R.
Shapley
A. E.
Kollmeier
J. A.
Zheng
Z.
Steidel
C. C.
Hainline
K. N.
2012
ApJ

745
33
Kunth
D.
Mas-Hesse
J. M.
Terlevich
E.
Terlevich
R.
Lequeux
J.
Fall
S. M.
1998
A&A

334
11
Kusakabe
H.
Shimasaku
K.
Nakajima
K.
Ouchi
M.
2015
ApJ

800
L29
Lake
E.
Zheng
Z.
Cen
R.
R.
Momose
R.
Ouchi
M.
2015
ApJ

806
46
Lawrence
A.
et al
2007
MNRAS

379
1599
Lilly
S. J.
et al
2009
ApJS

184
218
McCracken
H. J.
et al
2010
ApJ

708
202
McLinden
E. M.
et al
2011
ApJ

730
136
Martin
C. L.
Dijkstra
M.
Henry
A.
Soto
K. T.
Danforth
C. W.
Wong
J.
2015
ApJ

803
6
Matsuda
Y.
et al
2004
AJ

128
569
Matsuda
Y.
et al
2012
MNRAS

425
878
Matthee
J. J. A.
et al
2014
MNRAS

440
2375
Matthee
J.
Sobral
D.
Santos
S.
Röttgering
H.
Darvish
B.
Mobasher
B.
2015
MNRAS

451
400
Momose
R.
et al
2014
MNRAS

442
110
Nakajima
K.
et al
2012
ApJ

745
12
Nilsson
K. K.
Tapken
C.
Møller
P.
Freudling
W.
Fynbo
J. P. U.
Meisenheimer
K.
Laursen
P.
Östlin
G.
2009
A&A

498
13
Nilsson
K. K.
Östlin
G.
Møller
P.
Möller-Nilsson
O.
Tapken
C.
Freudling
W.
Fynbo
J. P. U.
2011
A&A

529
A9
Oesch
P. A.
et al
2015
ApJ

804
L30
Ono
Y.
Ouchi
M.
Shimasaku
K.
Dunlop
J.
Farrah
D.
McLure
R.
Okamura
S.
2010
ApJ

724
1524
Oteo
I.
Sobral
D.
Ivison
R. J.
Smail
I.
Best
P. N.
Cepa
J.
Pérez-García
A. M.
2015
MNRAS

452
2018
Ouchi
M.
et al
2004
ApJ

611
660
Ouchi
M.
et al
2008
ApJs

176
301
Ouchi
M.
et al
2010
ApJ

723
869
Ouchi
M.
et al
2013
ApJ

778
102
Overzier
R. A.
N. P. H.
Dijkstra
M.
Hatch
N. A.
Lehnert
M. D.
Villar-Martín
M.
Wilman
R. J.
Zirm
A. W.
2013
ApJ

771
89
Patrício
V.
et al
2016
MNRAS

456
4191
Reddy
N. A.
et al
2015
ApJ

806
259
Rivera-Thorsen
T. E.
et al
2015
ApJ

805
14
Rodighiero
G.
et al
2014
MNRAS

443
19
Rosdahl
J.
Blaizot
J.
2012
MNRAS

423
344
Rottgering
H. J. A.
R. W.
Miley
G. K.
van Ojik
R.
Wieringa
M. H.
1995
MNRAS

277
389
Sanders
R. L.
et al
2015
ApJ

799
138
Scarlata
C.
et al
2009
ApJ

704
L98
T.
et al
2014
ApJ

788
74
Skrutskie
M. F.
et al
2006
AJ

131
1163
Sobral
D.
et al
2009
MNRAS

398
L68
Sobral
D.
Best
P. N.
Matsuda
Y.
Smail
I.
Geach
J. E.
Cirasuolo
M.
2012
MNRAS

420
1926
Sobral
D.
Smail
I.
Best
P. N.
Geach
J. E.
Matsuda
Y.
Stott
J. P.
Cirasuolo
M.
Kurk
J.
2013
MNRAS

428
1128
Sobral
D.
Best
P. N.
Smail
I.
Mobasher
B.
Stott
J.
Nisbet
D.
2014
MNRAS

437
3516
Sobral
D.
et al
2015a
MNRAS

451
2303
Sobral
D.
Matthee
J.
Darvish
B.
Schaerer
D.
Mobasher
B.
Röttgering
H. J. A.
Santos
S.
Hemmati
S.
2015b
ApJ

808
139
Sobral
D.
Kohn
S. A.
Best
P. N.
Smail
I.
Harrison
C. M.
Stott
J.
Calhau
J.
Matthee
J.
2016
MNRAS

457
1739
Song
M.
et al
2014
ApJ

791
3
Steidel
C. C.
Erb
D. K.
Shapley
A. E.
Pettini
M.
Reddy
N.
Bogosavljević
M.
Rudie
G. C.
Rakic
O.
2010
ApJ

717
289
Steidel
C. C.
Bogosavljević
M.
Shapley
A. E.
Kollmeier
J. A.
Reddy
N. A.
Erb
D. K.
Pettini
M.
2011
ApJ

736
160
Swinbank
A. M.
Sobral
D.
Smail
I.
Geach
J. E.
Best
P. N.
McCarthy
I. G.
Crain
R. A.
Theuns
T.
2012
MNRAS

426
935
Swinbank
A. M.
et al
2015
MNRAS

449
1298
Taniguchi
Y.
et al
2015
ApJ

809
L7
Tasca
L. A. M.
et al
2015
A&A

581
A54
Taylor
M. B.
2005
Shopbell
P.
Britton
M.
Ebert
R.
Astronomical Society of the Pacific Conference Series Vol. 347, Astronomical Data Analysis Software and Systems XIV

Astron. Soc. Pac.
San Francisco
29
Trainor
R. F.
Steidel
C. C.
Strom
A. L.
Rudie
G. C.
2015
ApJ

809
89
Umehata
H.
et al
2015
ApJ

815
L8
Van Der Walt
S.
Colbert
S. C.
Varoquaux
G.
2011
preprint (arXiv:1102.1523)
Venemans
B. P.
et al
2007
A&A

461
823
Verhamme
A.
Schaerer
D.
Atek
H.
Tapken
C.
2008
A&A

491
89
Wardlow
J. L.
et al
2014
ApJ

787
9
Wisotzki
L.
et al
2015
preprint (arXiv:1509.05143)
Wold
I. G. B.
Barger
A. J.
Cowie
L. L.
2014
ApJ

783
119
Yang
H.
Malhotra
S.
Gronke
M.
J. E.
A.
Zheng
Z.
Dijkstra
M.
2015
preprint (arXiv:1506.02885)
Zheng
Z.
Cen
R.
Trac
H.
Miralda-Escudé
J.
2010
ApJ

716
574
Zitrin
A.
et al
2015
ApJ

810
L12