© ESO 2022

1. Introduction

An important question in galaxy evolution is how stars form in extremely metal poor environments. One way to address this is by observing star-forming galaxies at high redshifts. As already pointed out by Partridge & Peebles (1967), a tell-tale signature of such galaxies is their hydrogen Lyman α emission.

Young, massive, hot stars in star-forming regions produce hydrogen ionising radiation and when the hydrogen recombines, the most likely outcome of the recombination process is a Lyman α photon. Roughly ∼2/3 of the recombination events lead to the emission of such a Lyman α photon, assuming case B and a temperature of 10⁴ K (e.g., Dijkstra 2014). The higher the number of such massive, hot stars, the more ionising photons are produced, but stars that are bright in the Lyman continuum (LyC) are also the shortest lived, which means the starburst age, star formation rate, and stellar initial mass function (IMF) regulate the strength of the Lyman α radiation. A low metallicity in the gas forming the stars is typically expected to influence both the IMF, which is assumed to be universal except for metallicities of Z < 10⁻³ Z_⊙, and the number of ionising photons produced in such stars (Schaerer 2002; Raiter et al. 2010; Pallottini et al. 2015; Sobral et al. 2015).

Other effects to be considered for the production efficiency of ionising photons are the rotation of the stars (see e.g., Haemmerlé et al. 2017) and the evolution of binaries, which can strip stars of their envelopes resulting in hot, compact Helium stars which have a higher production of ionising photons (e.g., Stanway et al. 2016; Götberg et al. 2017), as well as the possible existence of super-massive stars with masses of M ∼ 10⁴ M_⊙ (Denissenkov & Hartwick 2014). Thus a young starburst with a high star-formation rate, top-heavy IMF, and low metallicity can result in the production of a high number of Lyman α photons. The ratio between the luminosity of the Lyman α line and the continuum luminosity density at λ_{Lymanα = 1215.67 Å} is defined as the equivalent width (EW, Schaerer 2003), which is usually corrected to the rest-frame equivalent width EW₀ = EW/(1 + z). This traces in principle the production rate of ionising photons per stellar mass. It is thus sensitive to the physical properties of a galaxy discussed above. Studying the distribution of EW₀ in Lyman α emitting galaxies, or Lyman Alpha Emitters (LAEs), therefore enables us to put statistical constraints on those properties.

It has often been quoted that there is an upper limit to the EW₀ of 240 Å (Charlot & Fall 1993; Malhotra & Rhoads 2002; Schaerer 2002; Laursen et al. 2013), although this value¹ is based on specific assumptions on the underlying stellar population (such as a constant SFR, an age of ∼10⁷ years, an IMF with a slope of 0.5, and an upper cut-off of the IMF of 80 M_⊙, Charlot & Fall 1993) and does not take into account radiation transfer effects through complex interstellar medium (ISM) kinematics and morphology, which can boost the EW₀ in addition. While a smaller rest-frame equivalent width can be explained by normal stellar populations or a lower escape fraction of Lyman α photons, a value exceeding this number implies unusual conditions in the star-forming regions of the galaxy on top of a high transmission through the ISM, circum galactic medium (CGM) and intergalactic medium (IGM) in the line of sight. Such unusual conditions can be sub-solar metallicities, high star-formation rates, young stellar ages, a top-heavy initial mass function, and a departure from case B (Schaerer 2003; Raiter et al. 2010; Maseda et al. 2020). While a young starburst even at solar metallicity can explain EW₀ ∼ 500 Å (Raiter et al. 2010), a low metallicity of Z/Z_⊙ ≲ 1/20, even for a normal IMF and older ages, can increase EW₀ by a factor of ≈70%. Thus selecting objects with large EW₀ is the most efficient way to find objects that are likely interesting to study when searching for example for young massive stars, low metallicity, or even population III stars.

Apart from photoionisation from stars and recombination, there are other channels that can produce Lyman α photons, even from outside the galaxy, which therefore do not correlate with stellar properties. In particular, gas falling in the potential well of galaxies, or shocks in the interstellar gas, can be cooled through de-excitation of collisionally excited hydrogen atoms, which produces Lyman α photons (a process that is not yet well understood, see reviews Dijkstra 2017; Faucher-Giguère et al. 2017). From numerical simulations, collisions are predicted to provide less than 10% of the intrinsic Lyman α emission produced in the galaxy, but could still represent 40% of the escaping radiation (e.g., Dayal et al. 2010; Mitchell et al. 2021). Observationally it is still difficult to determine the contribution of gravitational cooling of infalling gas to the production of Lyman α photons (e.g., Leclercq et al. 2017). Studies of Lyman α halos have pointed in a direction where this scenario is not favoured over the production in the star-forming regions of the galaxy or in satellite galaxies (Momose et al. 2016; Xue et al. 2017).

A third production channel is fluorescence, which is a photoionisation and recombination event where the ionising photon exciting the hydrogen atom was not produced inside the galaxy itself, but comes from the external ultraviolet (UV) background (e.g., Cantalupo et al. 2005; Kollmeier et al. 2010; Dijkstra 2017). This fluorescent emission always happens, but is predicted to be much fainter than the other production channels, except in the vicinity of a strong ionising source (Cantalupo et al. 2005). However, in so-called dark galaxies that are invisible in the UV but bright in Lyman α, it has been proposed that the main production mechanism can indeed be fluorescence, although stacking analysis shows that there is still ongoing star formation (e.g., Maseda et al. 2018). In both scenarios (collisional excitation and fluorescence), the strength of the Lyman α emission is not expected to relate to the strength of the UV continuum, and extreme EW₀ > 240 Å have been proposed as a way to identify dark galaxies, shining through Lyman α fluorescence triggered by a neighbouring quasar (Cantalupo et al. 2007, 2012; Marino et al. 2018).

On their way to the observer the Lyman α photons pass through the ISM, CGM and IGM. Since the photons are resonantly scattered by neutral hydrogen both in real and in frequency space, the strength and shape (e.g., Verhamme et al. 2006; Henry et al. 2015) of the Lyman α line will be affected by the neutral hydrogen column density (e.g., Shibuya et al. 2014a; Yang et al. 2017) as well as the dust content, the morphological properties, and kinematics of the gas (e.g., Atek et al. 2008; Trainor et al. 2016).

There are only a few effects that can potentially increase the observed EW₀, but most radiative transfer effects result in a decrease. For the former, there have been several studies concerned with the influence of radiative transfer effects on boosting the Lyman α equivalent width, for example through the preferential escape of Lyman α emission over the UV continuum (Neufeld 1991) or Lyman α beaming (Behrens et al. 2014; Behrens & Braun 2014; Verhamme et al. 2015). Numerical simulations have shown however, that the gain in EW₀ is marginal (e.g., Laursen et al. 2013; Gronke & Dijkstra 2014). It thus seems that large EW₀ values are more likely explained by the properties of the underlying stellar populations or through external effects like collisions and fluorescence (which would not influence the UV radiation), than through radiative transfer effects.

However, there are several reasons why we can expect the observed Lyman α EW₀ to be lowered. Since neutral hydrogen will cause the Lyman α photons to scatter, their path length is increased and thus also the probability to hit a dust grain and be destroyed or to get scattered out of the line of sight (see Dijkstra 2014 for a review), an effect that increases with increasing neutral hydrogen column density. The morphological and kinematic properties of the CGM can also play an important role in the radiative transfer of Lyman α photons and it has been shown that most Lyman α emitters have halos of extended Lyman α emission of ten times the size of the UV continuum (e.g., Hayashino et al. 2004; Steidel et al. 2011; Momose et al. 2014, 2016; Wisotzki et al. 2016, 2018; Leclercq et al. 2017, 2020). Given a high enough surface brightness sensitivity of the data, it is possible that such halos are even ubiquitous (see Kusakabe et al., in prep.).

The morphological structure of the galaxy itself and the angle at which we see it can also have an effect on the Lyman α line (e.g., Venemans et al. 2005; Gronwall et al. 2011; Bond et al. 2012; Jiang et al. 2013; Kobayashi et al. 2016; Paulino-Afonso et al. 2018). A clumpy morphology could be caused by a merger, multiple star-forming regions or satellite galaxies, each of which could boost the Lyman α emission. Paths driven by supernovae can allow Lyman α emission to escape in a preferential direction (Behrens et al. 2014), while the inclination of the galaxy can also lead to different EW₀ values (Verhamme et al. 2012; Behrens & Braun 2014), with face-on galaxies having potentially higher EW₀ (Laursen et al. 2009; Shibuya et al. 2014b).

The radiative transfer processes in the emitting galaxy can be imprinted on the shape of the Lyman α line, as a high neutral hydrogen column density in the ISM can broaden the Lyman α line or result in a blue bump (e.g., Verhamme et al. 2015). The neutral hydrogen in the IGM can then absorb any photons that are redshifted to the position of the Lyman α line wavelength, meaning any galactic emission to the blue of the Lyman α rest-frame wavelength, including the blue bump or the left side of the broadened line (e.g., Laursen et al. 2011; Garel et al. 2021), potentially changing its shape even after the photons left the galaxy.

Therefore, to get the Lyman α EW₀ directly after the galaxy, some studies assume a symmetric, usually single-peaked Lyman α line (after leaving the ISM and potentially CGM), the blue part of which would be attenuated by the IGM, which would lead to an apparent EW₀ of half the intrinsic value (e.g., Kashikawa et al. 2012; Zheng et al. 2014), where we define the intrinsic EW₀ as the rest-frame EW we would observe without IGM attenuation. The caveat with this idea is that first, the Lyman α line can be double peaked after passing through neutral hydrogen in the ISM or CGM (Hu et al. 2016; Matthee et al. 2018; Songaila et al. 2018; Meyer et al. 2021; Hayes et al. 2021) and contrary to the assumption that the IGM absorbs all emission to the blue of the intrinsic Lyman α wavelength, blue bumps at high redshifts have indeed been observed (e.g., Hayes et al. 2021). Second, it has been shown that the typical asymmetric line profile can be explained on the basis of the gas kinematics alone (Verhamme et al. 2006) and the Lyman α line is often shifted with respect to the systemic redshift (e.g., Kulas et al. 2012; Hashimoto et al. 2015; Verhamme et al. 2018; Muzahid et al. 2020), even before the emission reaches the IGM (e.g., in the low-redshift LAE analogues in the Lyman Alpha Reference Sample, see Östlin et al. 2014; Rivera-Thorsen et al. 2015, but also in Green Peas, see Henry et al. 2015). Although it is important to discuss the effects of the IGM on the Lyman α radiation when deriving the EW₀, it is not possible to correct for them for individual objects due to the redshift dependence and the high line of sight stochasticity of the opacity of the IGM (e.g., Inoue et al. 2014; Byrohl & Gronke 2020; Bassett et al. 2021).

In previous studies, different fractions of objects with large EW₀ have been determined among high-redshift LAEs. Typically, LAEs are found using the narrow-band technique (e.g., Hu et al. 1998; Malhotra & Rhoads 2004; Ouchi et al. 2008), which however needs spectroscopic follow-up observations to rule out the possibility of low redshift interlopers (e.g., Erb et al. 2014). There is no consensus yet on the fraction of high EW₀ values among high redshift LAEs, which is what we want to analyse in this paper.

Although purely narrow-band selected samples can produce a large number of LAE candidates, they need spectroscopic follow-up observations. In addition, they often lack the observational depth of deep Hubble Space Telescope (HST) data, since they get their large number of LAEs through their large surface area as they are restricted in the redshift range by the wavelength range given by the narrow-band. This makes it difficult to get good constraints on the UV continuum and thus on the EW₀ of Lyman α. What is more, Maseda et al. (2020) argue that surveys based on a pre-selection by UV continuum detections of LAEs in imaging are biased against finding high EW₀ objects.

Moreover, taking the full extent of Lyman α haloes into account when determining the EW₀ of Lyman α lines is important and often difficult for narrow-band surveys or even slit-spectroscopy. When measured to large enough radii and to lower surface brightness limits, many galaxies that otherwise would not have been classified as such, are actually LAEs (Steidel et al. 2011). The contribution of the Lyman α halo is typically around ∼65% (Leclercq et al. 2017), which means it often even dominates the Lyman α emission. For these reasons it is ideal to use integral field spectroscopy, to include the full Lyman α halo flux and correctly identify galaxies as LAEs.

Therefore we use data from the integral field spectrograph Multi Unit Spectroscopic Explorer (MUSE, Bacon et al. 2010) at the Very Large Telescope (VLT) in Chile, which is ideal for the detection of LAEs (e.g., Bacon et al. 2015; Wisotzki et al. 2016), combined with deep HST data to constrain the EW₀ as best as possible. MUSE has been used in previous studies of Lyman α EW₀ distributions, such as by Hashimoto et al. (2017a), who find a rather low fraction of high EW₀ > 200 Å LAEs (12 out of 417) using MUSE data of the Hubble Ultra Deep Field (HUDF, Beckwith et al. 2006; Bacon et al. 2015). This study mostly probes the low luminosity part of the LAE luminosity function, which is why we use data from both the MUSE-Wide (Herenz et al. 2017; Urrutia et al. 2019) and MUSE-Deep (Inami et al. 2017) surveys, containing a large sample of LAEs (around 2000) over a wide redshift range (3 ≲ z ≲ 6) and large survey area (100 pointings) of different observation depths.

We build on the work of Hashimoto et al. (2017a), who construct a Lyman α EW₀ distribution with the LAEs from the MUSE-Deep survey only and we add the MUSE-Wide detections to get better statistical constraints and probe a larger field-of-view and a larger range of Lyman α and UV luminosities.

This paper is structured as follows: In Sect. 2 we describe the data from the MUSE-Wide and MUSE-Deep surveys, as well as the ancillary deep HST broad-band data. In Sect. 3 we construct a sample of LAEs and describe the detection and classification of emission lines in MUSE data. We determine the UV continuum counterparts in the HST data and measure the UV continuum flux density as described in Sect. 4. In Sect. 5 we explain the line flux measurements as well as the fitting of the Lyman α lines. We combine the line fluxes measured from the MUSE data with the continuum flux densities measured from HST to obtain EW₀ in Sect. 6, where we show the EW₀ distribution and discuss the differences between the MUSE-Wide and the MUSE-Deep surveys. In Sect. 7 we show connections between EW₀ and other properties of the LAEs and in Sect. 8 we discuss the properties of the LAEs with the highest measured EW₀ and give a summary and conclusion in Sect. 9. Throughout this paper we use AB magnitudes, physical distances and assume flat ΛCDM cosmology with H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3 and Ω_Λ = 0.7.

2. Data

Measuring EW₀ can be dissected into two parts: determining the Lyman α line flux and the rest-frame UV continuum flux density. For both measurements we used different kinds of data, which made it possible to combine the power of integral field spectroscopy (for the line flux and other emission line properties) from MUSE (Bacon et al. 2010) with the depth of HST broad-band photometry (for the UV continuum). Therefore, special care had to be taken to properly combine the information gained from both data sets.

We used the spectroscopic data from MUSE to detect the LAEs and measure the line properties such as line flux, FWHM, asymmetry, and spectral peak separation (in case of a double peak in the spectrum, which is common for Lyman α lines). However, even for the MUSE-Deep data, the exposure time is not long enough to reliably measure the rest-frame UV continuum for all LAEs directly from the MUSE data, which is needed for the determination of the EW₀. Therefore, we used deep, broad-band data from HST at wavelengths longer than the redshifted Lyman α line so as not to contaminate the HST bands with the Lyman α emission itself.

2.1. Spectroscopic information from MUSE

MUSE is an integral field spectrograph and as such has two spatial and one spectral dimensions. It has a field of view (FoV) of one arcmin² and a spatial sampling of 0.2 arcsec. The spectral range covers 4750 Å–9350 Å which allows for the detection of Lyman α at 1215.67 Å in the redshift range of 2.9 < z < 6.7. Since this makes MUSE ideal to detect high-redshift LAEs, there are several surveys within the MUSE consortium which are aimed at studying their properties. In this paper we use data from the MUSE-Wide² (Herenz et al. 2017; Urrutia et al. 2019) and MUSE-Deep³ (Bacon et al. 2017; Inami et al. 2017) surveys. The MUSE-Wide survey, located in the CANDELS-Deep (Giavalisco et al. 2004; Koekemoer et al. 2011) as well as COSMOS (Scoville et al. 2007) fields, is specifically aimed at discovering a large number of bright LAEs around L^* of the LAE luminosity function (LF, see Herenz et al. 2019), while the MUSE-Deep survey samples more sub-L^* LAEs.

The two surveys, both taken at the ESO-VLT during guaranteed time observations (GTOs) of the MUSE consortium, complement each other, as they have different depths and sizes. The MUSE-Wide survey provides a shallower approach, with one hour exposure time for each field, but covers a large survey area with a total of 100 fields. Part of the data used in this paper is already publicly available, see the data release of the first 44 fields (Urrutia et al. 2019 which also contains a footprint of the MUSE-Wide survey). The pointings overlap by ≈4″, which reduces the total area covered to slightly less than 100 arcmin², covering a comoving volume of roughly ∼10⁷ Mpc³. The main part of MUSE-Wide is in the GOODS-S/CDFS and CANDELS-COSMOS areas with 60 fields. Eight additional fields are in the two HUDF parallel fields and 23 fields are in the COSMOS area. The nine HUDF fields of the MUSE-Deep survey area complete the 100 fields of MUSE-Wide (see Fig. 1 by Urrutia et al. 2019 for a footprint of the MUSE-Wide survey). The data reduction process for both MUSE-Wide and MUSE-Deep is described by Urrutia et al. (2019), using the data reduction pipeline by Weilbacher et al. (2020).

The detection limit for the Lyman α luminosity L_Lyα varies with wavelength due to the sensitivity curve of MUSE, the luminosity distance and the sky spectrum, as can be seen for example in Figs. 6 and 7 by Herenz et al. (2019), which show the selection function for MUSE-Wide LAEs in the first 24 fields. At a redshift of z ∼ 3 the 15% completeness limit is roughly log₁₀(L_Lyα[erg s⁻¹]) ≈ 42.2 and goes up to log₁₀(L_Lyα[erg s⁻¹]) ≈ 42.5 at a redshift of z ≈ 6.5 between the skylines. Herenz et al. (2019) find a characteristic Lyman α luminosity for their luminosity function of $\log L^{\ast }[\mathrm{erg}{\mathrm{s}}^{-1}]=42.2_{-0.16}^{+0.22}$.

The MUSE-Deep survey focuses on the HUDF, with the MUSE-Mosaic consisting of nine fields of ten hours exposure time each with a final contiguous area of 9.92 arcmin² (Bacon et al. 2017). The third and deepest tier in this construction is the ultra deep field (UDF) 10 (or MUSE-UDF10), located in the deepest part of the HUDF and reaching a total exposure time of 31 hours (see Bacon et al. 2017 for more information on the data reduction and Bacon et al. 2021 for the latest, deepest MUSE surveys).

We included the ten HUDF fields in this study (the nine mosaic fields and the deeper UDF10 field) to better understand the influence of the survey depth on the distribution of measured Lyman α EW₀ and to increase the sample size. Since the MUSE-Deep fields have a longer exposure time and thus go deeper in Lyman α luminosity and since the Lyman α luminosity function is steep, the number of objects in the 91 MUSE-Wide fields and in the ten MUSE-Deep fields including the UDF10 is similar, with 1051 objects from MUSE-Wide and 869 from MUSE-Deep. As the MUSE-Deep area has a different orientation and slightly overlaps with several other MUSE-Wide fields, objects that fall into both the MUSE-Wide and the MUSE-Deep parts of the survey were found by positional matching, and in case of duplicates only the MUSE-Deep information was used here.

From Fig. 5 of Drake et al. (2017a) the 15% completeness limit for the detection of Lyman α emission in MUSE-Deep is roughly log₁₀(L_Lyα[erg s⁻¹]) ≈ 41 at a redshift of z ≈ 3 and goes up to log₁₀(L_Lyα[erg s⁻¹]) ≈ 42 at z ≈ 6.5 (for the nine MUSE-Deep mosaic fields). The seeing conditions for MUSE-Wide were around $1\stackrel{″}{.}0$ for most of the MUSE-Wide fields (the full width at half maximum, FWHM, of the point spread function) and $\approx 0\stackrel{″}{.}6$ at 7750 Å for MUSE-Deep.

Combined, these three depths represent different Lyman α luminosities which allows for example to study the Lyman α luminosity function of LAEs, probing the bright end (with MUSE-Wide, see Herenz et al. 2019) as well as the faint end (with MUSE-Deep, see Drake et al. 2017a) but also the EW₀ distribution for UV-faint LAEs (Hashimoto et al. 2017a), which is a precursor to the present study. However, Hashimoto et al. (2017a) only consider LAEs in their sample that are detected with at least 2σ significance in at least two HST bands, while we also considered HST undetected objects.

2.2. Photometry from HST

We used the ACS F775W, ACS F814W, WFC3 F125W, and WFC3 F160W bands (Giavalisco et al. 2004; Koekemoer et al. 2011; Grogin et al. 2011) to measure the UV continuum flux density, the ACS F814W band to determine the UV continuum counterparts (and the ACS F775W band for the HUDF parallel fields where the ACS F814W was not available), and ACS F435W and ACS F606W to check for low-redshift interlopers⁴. This process is described further in Sect. 4.1.

3. Sample selection

In this section we describe the process of assembling a sample of LAEs from the MUSE-Wide and MUSE-Deep survey areas, including the detection of emission lines and their classification and report on a comparison with the catalogue by Inami et al. (2017).

3.1. Detection and classification

Integral field spectroscopy is ideal for detecting LAEs and we used the full two-dimensional spatial and combined spectral information of MUSE to detect emission line objects with LSDCat⁵ (Herenz & Wisotzki 2017), a tool using a matched-filtering approach to maximise the signal-to-noise (S/N) of the detections. For the first 24 fields we estimated the variance from 100 empty sky positions in the MUSE data directly (since the propagated noise from the pipeline underestimates the uncertainties due to the resampling when constructing the cube), but for the rest of the fields we used an updated version of the effective variance measurements (see Urrutia et al. 2019), taking into account the spectral resampling, which makes up a factor of 1.25 in the noise level, correcting the initial threshold for the first 24 fields of S/N = 8 to S/N = 6.4. For the following fields the threshold was lowered to five, which kept the rate of false positives at ≈5% when compared to the MUSE-Deep catalogue by Inami et al. (2017) (see explanation by Urrutia et al. 2019). This means that since the selection of LAEs in this paper was based on the detection of the Lyman α line in the MUSE spectra, we did not have a cut in EW₀ to define an LAE, but only in the S/N of the emission line. The full 100 MUSE-Wide fields with an updated data reduction and a consistent cut in S/N for emission line detection will be published by Urrutia et al. (in prep.).

After their detection, the emission lines were then grouped together by LSDCat into individual objects where the line with the highest S/N is referred to as the lead line. For consistency and to create a homogeneous sample, we used detections and line flux measurements from LSDCat for both the MUSE-Wide and MUSE-Deep fields (but using the full observational depth of the MUSE-Deep data). Therefore our sample of objects for MUSE-Deep is slightly different from that presented in the catalogue of Inami et al. (2017) and used for the EW₀ distribution study by Hashimoto et al. (2017a). An additional benefit of using LSDCat for measuring the Lyman α line fluxes is the inclusion of the full emission, even for objects with extended Lyman α halos (see the discussion on the influence of the halo size in Appendix A).

After detecting all emission line object candidates with LSDCat, the next step was to classify the objects, which was done with the help of the custom graphical user interface (GUI) QtClassify⁶, (Kerutt 2017, see also the Appendix of Herenz et al. 2017). QtClassify enables the user to load a full MUSE datacube, as well as a catalogue created by LSDCat and ancillary HST data if needed. For a detailed discussion of the procedure see Herenz et al. (2017) and Urrutia et al. (2019). For the nine MUSE-Deep fields and the UDF10 we used the same procedure as for the rest of the MUSE-Wide fields to create a homogeneous sample with the same detection and classification methods.

Even though MUSE provides spectra for each object, allowing to look for additional emission lines that could confirm or rule out the classification of an LAE, there is still the possibility of false classifications (see Urrutia et al. 2019 for a discussion), which can have an influence on the measured characteristic EW₀ values w₀ of the histograms of equivalent widths (see Sect. 6.1). To assess the reliability of the classification, we introduced a classification confidence between zero (the lowest confidence) and three (the highest confidence). This confidence was based on several aspects of the object, that each individual classifier used to base their decision on. If multiple lines matching a common redshift were detected in the spectrum, the confidence was three. If there was only one detected line but others were visible by eye (at a S/N below the detection threshold) using QtClassify, the confidence was set to two or three. If only one emission line was visible, the shape of the line was considered. Especially Lyman α lines often have a characteristic asymmetric and/or double peaked shape which is a clear indicator and resulted in a confidence of two to three. In the case of double peaks, we made sure that the peak separation and shape did not match with the O II doublet. If there was only one line, the S/N was low, and/or the line shape was not characteristic, a confidence of one was given. A confidence of zero was only given in extreme cases where no estimate of the redshift could be made. The process of classification was performed by two people who consolidated their results with a third person.

3.2. Comparison to the MUSE-Deep catalogue

For the MUSE-Deep part of our sample, a catalogue of objects using a different detection and classification method has been published by Inami et al. (2017) with a subsequent analysis of Lyman α EW₀ by Hashimoto et al. (2017a). In contrast to the MUSE-Wide approach this catalogue is not based solely on an emission line selection. Instead, several different methods were used to construct the catalogue, including an automated emission line detection software similar to LSDCat based also on a matched filtering approach (ORIGIN, Mary et al. 2020) and MUSELET (Piqueras et al. 2017) based on SExtractor (Bertin & Arnouts 1996), using narrow-band images created by collapsing five wavelength layers and subtracting the continuum (Inami et al. 2017). To classify the detected emission lines, the software MARZ (Hinton et al. 2016) was modified and the classification was done semi-automatically, aided by a team of human classifiers.

We matched the LAEs found with the MUSE-Wide approach to the LAEs in the Inami et al. (2017) catalogue for comparison, with a maximum allowed separation of $0\stackrel{″}{.}5$ and chose the closest counterpart within this area (given the redshifts in both catalogues matched as well). We chose this distance based on Herenz et al. (2017), where it was found that the 3σ positional uncertainty between the MUSE-Wide position and HST catalogues was below $0\stackrel{″}{.}5$. Of the 807 LAEs in the Inami et al. (2017) catalogue⁷, 124 could not be matched to detections with the MUSE-Wide method, while 214 objects were found in addition to the LAEs by Inami et al. (2017). The latter had on average a lower confidence level (with a mean of 1.7 compared to 2.1 for all LSDCat detected objects in MUSE-Deep), which means mostly weak lines were missed before. The 124 objects that were found in Inami et al. (2017) but could not be matched to our catalogue can be explained by the different catalogue creation methods. While Inami et al. (2017) not only searched for emission lines directly in the MUSE data, but also used the HST-based UVUDF catalogue (Rafelski et al. 2015) as priors, the MUSE-Wide method using LSDCat is a blind emission line search only⁸. Objects that are visible in HST but have emission lines below our detection threshold would therefore be missed. Conversely, LSDCat detected additional emission line objects that do not have an HST counterpart, which accounts for the 214 objects found with the MUSE-Wide method that are not in Inami et al. (2017). It should also be noted that we are only comparing LAEs here, which means that part of the discrepancy can be attributed to a difference in the classification, meaning LAEs that are missing in either catalogue might be present in the other but classified as something else.

3.3. Final sample of LAEs

We excluded three AGN from our sample, two of which (IDs 104014050 and 115003085) were already mentioned in Urrutia et al. (2019), one AGN (ID 214002011) is in the COSMOS field. For the MUSE-Deep part we excluded two objects (IDs 1841655 and 1381485), the first based on Hashimoto et al. (2017a), the second is a CIV emitter (both are also discussed in Bacon et al. 2021).

The total sample of LAEs used in this paper consists of 1920 LAEs. We excluded 35 objects that might be superpositions (two objects at different redshifts that overlap spatially), based either on a spectral energy distribution (SED) that does not match expectations for high redshift LAEs, which means no drop in the HST band to the blue of the Lyman α line (see Sect. 4.1) or on the presence of other emission lines in the MUSE data. In the latter case, we could not assign each emission line in the spectrum to one single redshift, which means that within the MUSE resolution, we were not able to disentangle the superposed objects. Since the reason for this is usually a low redshift interloper and to avoid any possible contamination, we excluded such objects, even if they were separated in the HST data. We thus got ∼11 LAEs per field in the MUSE-Wide survey (and therefore also roughly per square arcminute) and ∼96 LAEs per field for MUSE-Deep. To keep the methods consistent, we not only used the same detection method for MUSE-Wide and MUSE-Deep, but we also selected the UV continuum counterparts and measured the UV continuum flux density for MUSE-Deep in the same way as for the MUSE-Wide objects.

4. The UV continuum

As mentioned above, for measuring the EW₀ of the LAEs found in the MUSE data, the rest-frame UV continuum flux density was obtained from the deeper broad-band HST data instead of directly from the extracted MUSE spectra, since the latter were not deep enough for all objects. An additional advantage of the HST data is the higher resolution, allowing for a more detailed analysis of the UV continuum morphological properties. For this purpose we determined the UV continuum counterparts and fitted them with Galfit (Peng et al. 2002, 2010) where possible, a process which we describe in this section.

4.1. Identification of UV continuum counterparts

We determined the UV continuum counterparts of the LAEs in our sample in the HST data using the filter band ACS F814W for each object where possible (both for MUSE-Wide and MUSE-Deep, except for the eight HUDF parallel fields where we use the HST filter ACS F775W instead). We used the reddest available ACS filter band to have a high spatial resolution (twice better than the WFC 3 bands, see Table 1) and to make sure as many LAEs as possible would be detectable. Since the HST filter bands ACS F775W and ACS F814W overlap significantly (while ACS F814W is slightly deeper), using ACS F775W for the parallel fields was an adequate solution where ACS F814W was not available. We did not use any existing catalogues for the determination of the UV continuum counterparts as we wanted to be as unbiased as possible. Any signal within a radius of $0\stackrel{″}{.}5$ (measured from the maximum S/N in the Lyman α detection in MUSE) was taken into consideration as a counterpart (as this distance was found to be the 3σ positional difference by Herenz et al. 2017 when comparing the MUSE-Wide catalogue from LSDCat and the catalogue from Skelton et al. 2014). This criterion was used as a starting point, as we expected that not all Lyman α emission lines we find in the MUSE data have a UV continuum counterpart that is bright enough to be visible in the HST data. In case of more than one counterpart candidate within the $0\stackrel{″}{.}5$, additional HST filters were examined by at least two people to visually determine if the SED matches with what we would expect from LAEs, namely that there is little to no flux to the blue of the Lyman α line (as seen in the HST bands ACS F435W and ACS F606W), the spectrum declines towards the red, and there is an increase of flux in the band containing the Lyman α line.

Due to the IGM absorption, the emission will likely drop in HST filter bands to the blue side of the Lyman α line, although this depends on the assumed SED models. In cases where the possible counterparts were close and seemed to have the same positions and brightness in all HST bands considered, there is a possibility that it is really only one counterpart consisting of multiple components or clumps, which belong to the same object (see Fig. 1 for examples of different configurations). This could be caused by the clumpy nature of high redshift galaxies, due to distinct star formation areas or satellite galaxies. In this case, all components (within $0\stackrel{″}{.}5$) were considered to belong to the same Lyman α emission (as e.g., in the left panel of Fig. 1). When matching the Lyman α positions to the UV continuum counterparts, there is often a slight shift between the two positions (see also Claeyssens et al. in prep. for a study using lensed LAEs), as the Lyman α emission might escape preferentially through outflows or holes in the ISM of the galaxy (e.g., Shibuya et al. 2014b).

Fig. 1. Examples of LAEs with different kinds of UV continuum counterparts, shown in 2″ × 2″ cutouts of the filter band ACS F814W. The white crosses are centred on the UV continuum position, the white circles are centred on the Lyman α position from LSDCat and have a diameter of $0\stackrel{″}{.}5$. Left image: object with multiple components or clumps in the counterpart. Middle image: object without a significant detection in the filter band ACS F814W (the centre of the cross is set on a noise peak/bright pixel). Right image: counterpart with only one component, but slightly offset with respect to the Lyman α position.

We show the redshift distribution of our sample in Fig. 2 (we use redshifts corrected using the equations in Verhamme et al. 2018 based on the information of the FWHM and peak separation of the line, see Sect. 5.2). The constant decrease in the number of LAEs with high redshifts is influenced by the increasing luminosity distance and the resulting detection of more luminous objects, which are less numerous (see Herenz et al. 2019 for a discussion on the selection function and Lyman α luminosity function). At redshifts z > 5 the number of objects with and without a UV continuum counterpart is almost the same, while at lower redshifts most objects can be seen in the HST data. This is caused by the decline in the intrinsic UV continuum brightness at higher redshifts as well as cosmological surface brightness dimming, making it harder to detect the counterparts, commonly referred to as the Malmquist bias (Malmquist 1922, 1925), which also contributes to the decline in LAE numbers at higher redshifts.

Fig. 2. Redshift histogram of the full sample of LAEs, split into groups according to their UV continuum counterparts in HST, in bins of Δz = 0.25. The dark blue histogram shows objects with UV counterparts which are not resolved into multiple components. The turquoise histogram shows objects without a visible UV counterpart in the available HST imaging and the light green histogram shows objects with UV counterparts consisting of more than one component. In these cases, the HST counterpart could not be fitted with a single Sérsic profile, as it was clumpy or asymmetric.

4.2. Fitting the UV continuum using Galfit

After assigning UV continuum counterparts to all LAEs where possible, we fitted them with Galfit, a fitting algorithm created for HST images that is fast and efficient at fitting multiple components simultaneously (Peng et al. 2002, 2010). From these fits we gained information not only on the magnitudes, but also on the sizes, number of components, and axis ratio of the UV continuum counterparts where they are resolved in HST. Another advantage of using Galfit over simple fixed apertures is the possibility to model (and thus subtract) neighbouring objects as well, to obtain more reliable magnitude measurements. As the LAEs we want to fit are at high redshifts and thus small even in the HST data, we used Galfit to fit them with a simple Sérsic (Sersic 1968) profile (per component):

$$\begin{array}{c}\hfill \log \left(\frac{I(R)}{I_e}\right)=-b_n[{\left(\frac{R}{R_e}\right)}^{\frac{1}{n}}-1].\end{array}$$(1)

With I(R) the surface brightness, I_e the surface brightness at R_e (the effective radius containing half of the luminosity), b_n ≈ 2n − 0.331 (Peng et al. 2002) making sure that half of the total flux will be within one effective radius, and n the Sérsic (power-law) index. Depending on the Sérsic index n, the profile can become a Gaussian for n = 0.5, an exponential profile at n = 1, and a de Vaucouleurs profile (de Vaucouleurs 1948) with n = 4.

To facilitate the fitting, we used a wrapper with a GUI that automatically generates input files for Galfit, first based on initial priors⁹, then on previous runs, and stores the output data in a coherent way (see Fig. 3 for a screenshot of the wrapper). We iterated the fitting using the results form the previous run as new initial priors until the parameters were stable, in cases where the fitting was difficult. The procedure was to fit the morphological properties of the objects first in the ACS F814W filter (or the ACS F775W for the HUDF parallel fields), which is the deepest in the MUSE-Wide area, and then use them as fixed values for the other filter bands (ACS F606W, ACS F775W, WFC3 F125W and WFC3 F160W) with the magnitude as a free parameter. That way fainter flux at larger wavelengths could be captured within the same area as in the deepest HST filter. Another advantage of fixing parameter values was that even if the object was faint and its morphological parameters could not be fit well with Galfit, we could capture the magnitude by reducing the number of free parameters, facilitating the fitting. For faint objects, the axis ratio and Sérsic index n could often not be constrained in the fit. In this case we fixed n = 1 and b/a = 1, that is we fitted only a circular exponential profile, so that at least the effect radius and the objects’ magnitude are recovered. In some cases the objects were still too small to be resolved in HST and we used a fixed aperture to capture the magnitudes. Of the 1284 objects that have an HST counterpart, 538 have a reliable measurement of the effective radius.

Fig. 3. Screenshot of the Galfit wrapper used for fitting and measuring the UV continuum in different bands. The left image shows a cutout of the data in the HST filter band ACS F814W, the middle image shows the (in this case) two models for the object and the right image shows the residuals. Clicking in the right panel will allow the user to add a Galfit model for an additional component, which creates a green circle to mark the place. Below the panels are morphological measurements as well as magnitude and position, each of which can be fixed if needed. In the example, the object was fitted with two Sérsic models, each of which can be selected by clicking on the small green circles. The check boxes on the right allow the user to asses whether the object is clumpy and if there are multiple components in the counterpart or no counterpart at all, as well whether to include the model in the object.

In cases where the object was very faint or maybe even not visible, we also used a fixed radius of the size of the point spread function (PSF) for the Galfit model. If the flux in this aperture was not above 1σ in any of the HST bands, we considered the object undetected in HST.

For the HST PSF we used a Moffat function (Moffat 1969) measured using Galfit, from stars in the corresponding HST filter bands. With the help of the wrapper, nearby objects that might influence the fit of the LAE can be fit as well, so they do not artificially increase the measured magnitude of the LAE. In the same way, multiple components can be fit and included in the total flux of the UV continuum counterpart. In these cases, we added the continuum flux density of all components for the EW₀ estimation, assuming that the Lyman α emission is produced from the entirety of the components. To exclude the possibility of multiple components being superpositions of different objects, we checked for additional emission lines in the MUSE spectra and for different SEDs in the separate components of the counterparts in the HST data (see Sect. 4.1 above). For the morphological properties of the LAEs, that is the effective radius and axis ratio, we used the parameter uncertainty estimates from the Galfit models, while we used random apertures for the errors on the magnitudes (see Sect. 4.4 below).

Using this method of carefully checking and measuring the UV continuum counterpart properties of each individual object we can go beyond existing catalogues (e.g., Guo et al. 2013; Skelton et al. 2014; Rafelski et al. 2015) and include even very faint sources that had been missed before, but that we know to be real through the Lyman α line we found with MUSE. If we compare objects for which we did find an HST counterpart to the catalogues by Guo et al. (2013) and Skelton et al. (2014), we find that 42% of our LAEs have a match in the Guo et al. (2013) catalogue and 59% have a match in the Skelton et al. (2014) catalogue. The discrepancy can be explained by the fact that these catalogues are based on near-IR detections, while our LAEs are most likely UV-dominated.

4.3. Continuum flux densities

For the EW₀, we need the continuum flux density at the position of the Lyman α line (for objects where a counterpart is detected in the HST data). In principle the flux density at λ_Lyα could be obtained directly from the HST band containing Lyα. However, this requires a correction for the contribution of the Lyman α flux to the band, which introduces uncertainties that are hard to quantify. This is because of the halo component that is often undetected at the depth of HST. Moreover, the necessary correction for IGM absorption to the blue side of Lyman α introduces an additional complication, since we do not know how much of the measured flux in the HST filter band can be attributed to the UV continuum redwards of the line. Therefore we decided not to attempt a correct for the Lyman α line flux or the IGM absorption and instead we used the flux density from the HST filter band to the red of the line that has zero throughput at the Lyman α wavelength.

This means for objects below a redshift of ∼4.7 we used the HST bands ACS F775W, ACS F814W, WFC3 F125W and WFC3 F160W. For objects above this redshift we used the HST bands WFC3 F125W and WFC3 F160W to measure the magnitudes which are then used to estimate the flux at the Lyman α position.

To extrapolate the measured flux at the effective wavelength of the HST filter to the Lyman α line position, we need to know the UV continuum slope (β) of the spectrum, which correlates the flux density f at a certain wavelength λ to the wavelength via f_λ ∝ λ^β. If we know the flux density f at two or more different wavelengths λ₁ and λ₂, we can derive the UV continuum slope, assuming the continuum is a power law:

$$\begin{array}{c}\hfill \beta =\frac{{log}_{10}\left(\frac{f_{{\lambda }_1}}{f_{{\lambda }_2}}\right)}{{log}_{10}\left(\frac{{\lambda }_1}{{\lambda }_2}\right)}.\end{array}$$(2)

If we have detections in two or more HST filter bands, we can fit a linear relation to the logarithm of the measured flux density values. For this we used HST bands to the red side of the emission line (as mentioned above) and fitted a simple linear relation to the continuum flux density measurements. For objects that have a detection in only one HST band or even no counterpart at all, it is not possible to measure the β parameter. For this reason, and since the measured β values scatter significantly and have large errors especially for faint objects, we used our median value of β = −1.97 of the entire sample as a fixed value for each individual object.

The caveat with this approach is that LAEs with high EW₀ possibly have different properties from LAEs with lower EW₀, which means their β values could be systematically different from each other. The same could apply for LAEs that are undetected in the UV continuum for which it is not possible to measure the continuum slope at all from individual objects. If such objects had a bluer continuum slope due to less dust or younger ages, that would mean that their real EW₀ would be lower than what we assume here (see Maseda et al. 2020). Therefore the lower limits we are quoting here for EW₀ of HST undetected objects are only lower limits assuming a fixed β = −1.97. Ideally, for objects that are individually undetected in the UV continuum, one could stack the HST photometry to obtain β values from the stacks (see Maseda et al. in prep.).

In the study on EW₀ distributions using data from MUSE-Deep alone, Hashimoto et al. (2017a) measured the UV continuum slope from three (or two) adjacent HST bands where available. They show that using a fixed value for β has little impact on the derived characteristic EW₀ values w₀, but caution that at lower (higher) redshifts, the EW₀ values can be overestimated (underestimated) when using a fixed β value of β = −2 due to the redshift evolution of the UV continuum slope.

Most literature studies find β values around β ≤ −2 (Bouwens et al. 2009; Castellano et al. 2012) for Lyman Break Galaxies (LBGs), but Karman et al. (2017) and Santos et al. (2020) find even steeper slopes for faint LAEs at similar redshifts. A larger UV continuum slope indicates a redder spectrum, corresponding to dust absorption and metallicity in the galaxy. The median of the UV continuum slope derived from multiple HST bands β = −1.97 used here thus matches well with results found in the literature. To understand the influence of different β parameters on the measured characteristic EW₀ values of the EW₀ histograms, we compare the median value to the upper β = −1.57 and lower β = −2.29 quartiles of the distribution of measured β values (see Sect. 6.1).

4.4. Limits and magnitude errors

Any object not detected in a particular HST filter above a 1σ detection significance was assigned the limiting flux density (assuming a point source) in that HST band as the continuum flux density (see Table 2 for median values of the limiting magnitudes in each band) and the HST filter was excluded from the estimations of EW₀ and the UV continuum slope. Objects with no significant detection in any HST filter band have lower limit EW₀ which we determined using the limiting flux in the deepest band (ACS F775W for the parallel fields and ACS F814W for all other bands) in combination with the median UV slope β = −1.97. In this paper, in plots showing the full sample including objects with lower limits for EW₀ (which is always the case except if explicitly mentioned in the figure caption), those limiting values are shown as if they were detections at their respective limits.

Since we used different fields that have different HST depths and since the depth can also vary over the field, we measured the limiting flux density for each object using 100 random apertures of $0\stackrel{″}{.}5$ close (within a radius of 2″) to the Lyman α position of the LAE in question. To avoid neighbouring objects, we excluded apertures with a measured flux density above the standard deviation of the 100 apertures as well as apertures that are outside the field of view of the specific MUSE pointing. In order to be consistent with the error measurements of the magnitudes, we used the same method for the 1σ errors on the flux densities of objects that do have measureable counterparts in the HST photometry.

4.5. Results: objects with no HST counterpart

In total, of our sample of 1920 LAEs, 636 (around 33% of the full sample) have either no visible counterpart in the UV continuum or only a faint one that was not distinguishable from noise (see Bacon et al. 2017; Marino et al. 2018; Maseda et al. 2018, 2020 for studies of dark galaxies not visible in the rest-frame UV continuum). The 869 MUSE-Deep LAEs have fractionally more UV continuum non-detections (44%) than the 1051 MUSE-Wide LAEs (24%). This is mainly caused by the increase in emission line depth in the MUSE data of MUSE-Deep, making it easier to detect objects in Lyman α that are too faint to see in the HST photometry (which does not have a big difference in depth in the different fields, see Table 2). This can be seen in Fig. 4, where the Lyman α luminosity histograms of MUSE-Wide and MUSE-Deep objects with and without counterparts are shown. The MUSE-Deep objects occupy the fainter range of Lyman α luminosity values and the objects without UV continuum counterparts have a slight shift to lower line luminosities as well (see Sect. 6.2).

Fig. 4. Histogram showing the logarithmic Lyman α luminosity of the LAEs in the MUSE-Wide (light green, labelled MW) and MUSE-Deep surveys (dark green, labelled MD). The samples are divided into objects with UV continuum counterparts and without counterparts. The number of objects per sample are stacked, which means every histogram bar gives the total number of objects in that luminosity bin.

4.6. Discussion: morphological properties of the UV continuum

Most objects (1118 or 58%) have a counterpart with only one component visible in the HST data, but 166 (≈8.6%) have a counterpart consisting of multiple components (see Fig. 2 and left panel of Fig. 1). LAEs consisting of multiple components could be intrinsically clumpy or irregular, which could have several causes: It could be that (1) a merger is underway, which could potentially boost the star formation and thus the Lyman α EW₀, (2) multiple star-forming regions give an irregular shape and outshine the rest of the galaxy or (3) multiple objects or satellite galaxies at close range that all contribute to the Lyman α emission (e.g., Venemans et al. 2005; Gronwall et al. 2011; Bond et al. 2012; Jiang et al. 2013; Kobayashi et al. 2016; Paulino-Afonso et al. 2018). Galaxies at high redshifts do not necessarily follow the same morphological patterns we see in the nearby universe and Gronwall et al. (2011) show that while LAEs are heterogeneous in their rest-frame UV morphological properties, they have more in common with high-z LBGs than with local galaxies and tend to be on average rather compact. However, the smaller sizes of LAEs could be the result of a possible redshift evolution, which is found by Shibuya et al. (2019), Ferguson et al. (2004), but for example not by Paulino-Afonso et al. (2018). Even though most of our LAEs are not resolved in HST, which could be caused by their compact nature, it is interesting to see a significant fraction of LAEs at redshifts 3 < z < 6 being clumpy or possibly merging.

5. The Lyman α line

Another ingredient to measuring the EW₀ is the line flux of the Lyman α line (Sect. 5.1). Its spectral properties can give us information on possible radiative transfer processes and the neutral hydrogen column density and we find a large variety of different Lyman α line shapes, as can be seen in Fig. 5, which we fitted with an asymmetric Gaussian (Sect. 5.2).

Fig. 5. Overview of the Lyman α emission lines of the 35 LAEs with the highest measured EW0, sorted by EW0 in descending order (left to right, top to bottom). We note that the y-axis shows the flux normalised to the maximum of each peak for better comparison. The vertical dashed lines indicate the positions of the Lyman α lines and the blue bump if there is one. The grey shaded area shows the standard deviation from the MUSE spectra. The shown objects are taken both from the MUSE-Wide and -Deep surveys and have been continuum subtracted.

5.1. Line fluxes

The line flux was measured directly from the MUSE data cubes using LSDCat (Herenz & Wisotzki 2017) in spatial apertures of three Kron radii (Kron 1980, see Herenz et al. 2017), which gives the first-order moment of the light distribution (Bond et al. 2009). The spectral range for the line flux measurement of LSDCat includes all voxels (volume pixels) above the analysis threshold (see Herenz et al. 2017, especially Fig. 8). This way we used the full three dimensional information from the MUSE cube to include the entire Lyman α line flux. In Appendix A we show a comparison between our measurements for line fluxes and measurements from Leclercq et al. (2017) including a modelling of the extended Lyman α halo. Since both values match well, we conclude that our line flux measurements capture well the flux in the extended Lyman α halo (except for extreme cases).

This extended emission is often not taken into account when measuring Lyman α EW₀ if the line flux is measured from slit spectroscopy or narrow-band images where the object is assumed to have the same size in Lyman α as in the complementary broad-band images. Since 40 − 90% of the line flux can be found in the Lyman α halo when comparing the Lyman α emission to the extent of the UV continuum (Wisotzki et al. 2016), a large portion of the line flux could be omitted. In cases where the Lyman α line is double peaked (see Sect. 5.2 below), LSDCat sometimes does not include the full flux of the blue bump, which is why we corrected the line flux measurements for the fraction missed in this way. For this correction we fitted a linear combination of two asymmetric Gaussian functions to one-dimensional spectra (see Sect. 5.2 below), took into account the wavelength window that LSDCat used for the line flux, measured the ratio of the part of the line that was missed to the part that was included, and corrected the line flux accordingly.

Another potential correction of the Lyman α line flux is the absorption in the IGM (e.g., Stark et al. 2011; Laursen et al. 2011; Caruana et al. 2014; Kusakabe et al. 2020; Hayes et al. 2021). While many studies on Lyman α EW₀ assume a generic flux correction, the line shape can be influenced by radiative transfer processes and gas kinematics even before absorption in the IGM. In addition to this, the stochasticity and redshift-dependence of the IGM absorption (see e.g., Thomas et al. 2017, 2020 who show the large dispersion in the IGM transmission along different lines of sight) make it impossible to correct at an object by object basis, also given the fact that LAEs could reside in ionised bubbles (e.g., Roberts-Borsani et al. 2016; Castellano et al. 2016, 2018; Stark et al. 2017; Mason et al. 2018; Tilvi et al. 2020; Jung et al. 2020; Endsley et al. 2021) and we do not know the systemic redshifts. Therefore we did not correct for IGM absorption, so as not to overestimate the Lyman α line flux and thus the EW₀. In extreme cases, assuming the Lyman α line is symmetric around the systemic redshift (which is unlikely though due to radiative transfer processes) and the IGM absorption in the sightline is 100%, we could thus underestimate the real Lyman α flux by half.

5.2. Lyman α line shape properties

In this section we quantify the line shape properties, which are the asymmetry, FWHM (corrected for the line spread function, LSF, of MUSE), double peak fraction, and peak separation from one-dimensional spectra extracted from the MUSE datacubes. These spectra were obtained by weighted summation in each spectral layer. As weighting function we used the wavelength dependent Moffat (Moffat 1969) profile that describes the PSF of our observations (see Urrutia et al. 2019). This was done to get the highest possible S/N in the emission line for our spectral fits, by weighing down the more noisy outer parts of the emission. This way we did not retain the information of the total line flux in these spectra, which is why we used the line flux measurements directly from LSDCat. We also did not model any specific spatially extended shape of the Lyman α emission for these one dimensional spectra that we used for fitting the Lyman α line. This assumes that there are no spatial variations in the line shape properties and the line has the same properties in the halo as in the central part of the LAE. This simplification is sufficient for the purpose of this paper, but it should be noted that a spatial variation of Lyman α line properties of our high-redshift LAEs is possible (see e.g., Erb et al. 2018 for an example at z = 2.3, Claeyssens et al. 2019 and Leclercq et al. 2020 for Lyman α halos in the MUSE Lensing Clusters survey, Richard et al. 2021, and the MUSE-Deep survey).

To determine which Lyman α lines have a blue bump, we used a visual inspection tool that also fits the spectra with an asymmetric equation (explained below, Eq. (3)). We opted for a manual detection of blue bumps for this work to find even unusual multiple peaks, but an automatic determination of double peaks is a possibility as well (see Vitte et al., in prep.). As mentioned above, we used the fits to the two parts of the Lyman α line to correct the line fluxes to include the blue bumps, since the line flux is taken from LSDCat in a certain wavelength window, which not always includes the full blue bump, especially for high peak separations.

In order to gauge the reliability of the visual inspection regarding the presence of blue bumps, we used a Monte-Carlo type approach by creating 1000 randomised spectra with varying S/N values (of the red, main peak compared to the noise). The basis for this were ten spectra with intrinsically high S/N ratios (all with a S/N > 20 for the main, red peak), five of them with clear blue bumps and five without blue bumps, taken from the first 24 fields of the MUSE-Wide survey. Their peak separations range from 308 ± 39 km s⁻¹ to 477 ± 89 km s⁻¹ and their blue bump to total line flux ratios range from 9% to 29%. These spectra were then artificially degraded to lower S/N values between 0.5 and 10 of the main line and again analysed in the same way as the real spectra, which means determining visually whether an object has a double peak or not. However, since we know from the original spectra which objects should have double peaks, we can now analyse up to which S/N value we are able to accurately recover the double peaks.

At the S/N values used as the detection limit in MUSE-Wide (6.4 and 5 respectively for the first 24 fields and the rest of the fields), we already reach an accuracy of over 80%, which means in 80% of cases we could correctly determine whether the spectrum has a blue bump or not (where the accuracy is the sum of determined true positives and true negatives divided by the sum of actual positives and negatives). This also depends on the peak ratio, since a low blue bump to main peak ratio in combination with a low S/N of the main peak will result in a blue bump that is harder to detect than for the same S/N with a higher peak ratio. Another aspect is the peak separation of the chosen objects. We chose objects that had obvious double peaks, which lead to the rather narrow range of peak separations. For closer peaks, the accuracy of the classification is expected to be lower.

We fitted the Lyman α line using the asymmetric Gaussian function described by Shibuya et al. (2014a):

$$\begin{array}{c}\hfill f(\lambda )=Aexp(-\frac{{(\lambda -{\lambda }_0^{\mathrm{asym}})}^2}{2{\sigma }_{\mathrm{asym}}^2})+f_0.\end{array}$$(3)

Here, f₀ is the continuum level, A is the amplitude, and ${\lambda }_0^{\mathrm{asym}}$ is the peak wavelength of the line. For the latter, we took the LSDCat measurements as a first guess. The asymmetric dispersion is σ_asym, consisting of ${\sigma }_{\text{asym}}=a_{\text{asym}}(\lambda -{\lambda }_0^{\text{asym}})+d$. Here, d is the typical width of the emission line and a_asym is the asymmetry parameter. A positive asymmetry value suggests a line with a red wing, which is the case for most of the red (main) Lyman α lines, a negative value means the line has a blue wing. As described above, a Lyman α line can also be double peaked. If that is the case and the blue bump was fit individually, the asymmetry only refers to the red peak. If there is a close but unresolved blue bump, the measured asymmetry will be smaller than if the double peak could have been resolved and fit individually or there could even be a blue wing. After fitting the line, we derived the FWHM value from the asymmetric Gaussian fit, which is given by (see also Claeyssens et al. 2019):

$$\begin{array}{c}\hfill \mathit{FWHM}=\frac{2\sqrt{2\ln 2}d}{1-2\ln 2a_{\mathrm{asym}}^2}.\end{array}$$(4)

We corrected the FWHM of the Lyman α line for the spectral line spread function (LSF) of MUSE, which can be approximated by a Gaussian. The FWHM of the LSF is wavelength dependent and we used the value for the MUSE-Deep Mosaic fields (with ten hours exposure time) given by Bacon et al. (2017), which follows F_mosaic(λ [Å]) = 5.835 × 10⁻⁸λ² − 9.080 × 10⁻⁴λ + 5.983. This correction assumes that the Lyman α line can be approximated by a Gaussian as well, which is not always the case. Therefore the LSF correction itself is also just an approximation (for a detailed discussion of this problem see Childs & Stanway 2018). It should be kept in mind as well, that the measure of the asymmetry becomes unreliable for narrow lines which are dominated by the LSF.

Verhamme et al. (2015, 2017) have predicted and shown that the line shape properties are connected with the Lyman α escape fraction and EW₀ of the LAEs. There is a correlation between the peak separation and the shift of the red peak with respect to the systemic redshift as well as between the FWHM of the line and the shift of the red peak (Verhamme et al. 2018). Thus, we used both the peak separation and the FWHM to more accurately estimate the systemic redshift of the LAEs (according to equations one and two in Verhamme et al. 2018). We use these corrected redshifts throughout the paper.

6. Equivalent widths

To investigate the strength of the Lyman α lines and later compare it to other properties, we measure EW₀. With the measurements of the Lyman α line flux $F_{\mathrm{Ly}\alpha }^{\mathrm{line}}$ and UV continuum flux density $f_{\mathrm{Ly}\alpha }^{\mathrm{cont}}$ at the wavelength of the Lyman α line described above, we can now determine the Lyman α EW₀ as the fraction between the two:

$$\begin{array}{c}\hfill \mathrm{EW}={\displaystyle \underset{{\lambda }_0}{\overset{{\lambda }_1}{\int }}\frac{f_{\mathrm{Ly}\alpha }^{\mathrm{line}}-f_{\mathrm{Ly}\alpha }^{\mathrm{cont}}}{f_{\mathrm{Ly}\alpha }^{\mathrm{cont}}}\mathrm{d}\lambda \approx \frac{F_{\mathrm{Ly}\alpha }^{\mathrm{line}}}{f_{\mathrm{Ly}\alpha }^{\mathrm{cont}}}.}\end{array}$$(5)

Here, λ₀ and λ₁ define the range of integration, which is the width of the Lyman α line, while the flux density in the line is $f_{\mathrm{Ly}\alpha }^{\mathrm{line}}$. The approximation used here on the right side of the equation is valid given $f_{\mathrm{Ly}\alpha }^{\mathrm{cont}}<<f_{\mathrm{Ly}\alpha }^{\mathrm{line}}$.

The rest-frame EW is then given as EW₀  =  EW / (1 + z). As explained above, we did not correct the measured EW₀ for potential IGM absorption. For all plots showing EW₀ measurements (except the histograms) we include both the error on the Lyman α line flux as well as on the UV continuum measurements in the error bars.

6.1. Results: histograms of equivalent widths

Histograms of EW₀ are shown in Fig. 6 with exponential fits to determine the scaling factor. Objects where only a lower limit for the EW₀ could be measured are shown in the histogram as the green bars. The best-fit values (using a least-squares fit) for the exponential fit N = N₀ exp(−EW₀/w₀) are w₀ = 75.8 ± 1.9 Å and w₀ = 95.5 ± 3.5 Å for the characteristic EW₀, and N₀ = 804.2 ± 20.3 Å and N₀ = 901.4 ± 5.6 Å (full sample and for the sample excluding lower limits, respectively). Results in the literature on w₀ vary but generally fit well with our measurements. Recently, Jung et al. (2018) found w₀ ∼ 60 − 100 Å for a sample of LAEs with redshifts 0.3 < z < 6. In a previous study on LAEs in the MUSE-Deep survey, Hashimoto et al. (2017a) divided their sample in three redshift bins, finding w₀ = 113 ± 14 Å, w₀ = 68 ± 13 Å and w₀ = 134 ± 66 Å for redshifts z ∼ 3.6, z ∼ 4.9 and z ∼ 6.0, respectively.

Fig. 6. Histograms of the Lyman α EW0 of the full sample, including both MUSE-Wide and MUSE-Deep objects. The blue histogram shows the objects with a secure measurement (excluding lower limits), the green bars on top show the objects with a limiting EW0. The exponential fits N = N0 exp(−EW0/w0) are shown as the dashed lines, blue for the objects with secure measurements and light green for the full sample including lower limits. The bin width in EW0 is 50 Å. For the measurement of the characteristic EW0w0 (given below the legend and in Table 4) the smallest EW0 bin was excluded as it is likely incomplete (as Lyman α lines with small or zero line fluxes are hard to detect and preferentially result in small EW0). The shaded areas around the dashed lines indicate the exponential fits for the EW0 distributions using UV continuum slopes of β = −1.57 and −2.29, the first and last quartiles of the distribution of measured β values.

Table 3 gives an overview of the number of objects in the full sample as well as in the MUSE-Wide and MUSE-Deep samples that have EW₀ > 100 Å and EW₀ > 240 Å. The strength of the MUSE-Wide survey in detecting extreme LAEs is evident when looking at the statistics: including the lower limits for EW₀ (for objects that are undetected in the HST data), around 20% of LAEs in MUSE-Wide have EW₀ higher than predicted by stellar population models (which set the limit at 240 Å, assuming solar metallicity, a constant SFR, and an upper cut-off of the IMF of 80 M_⊙, Charlot & Fall 1993), while only 11% of objects in MUSE-Deep have EW₀ > 240 Å.

To estimate the influence of the choice of β on the EW₀ histograms and the characteristic EW₀, we used the first and last quartile values of our measured UV continuum slopes in addition to the median value of β = −1.97 and show the result as the shaded range in Fig. 6 and in Table 4. This highlights that, within a reasonable range, the choice of UV continuum slope does not affect the measured characteristic EW₀ much.

By using the largest, homogeneous dataset of MUSE-identified LAEs, combined with the deepest HST broad-band data, we are able to make accurate measurements of EW₀ for ∼2000 galaxies, making this sample over an order of magnitude larger than previous studies. Our measurements establish the existence of high EW₀ objects with EW₀ > 240 Å and even several 10² Å, and their occurrence rate is not low, so the physical conditions allowing for the production of such a huge number of Lyman α photons per UV magnitude seem rather common in 3 < z < 6 star-forming galaxies. Other studies often correct for the Lyman α absorption of the IGM, but as explained in Sect. 5.1, we did not correct for the IGM. This means that some of our measured EW₀ could intrinsically be even larger (if we assume that the intrinsic EW₀ is what comes out of the galaxy without IGM attenuation).

Gawiser et al. (2006) caution that non-detections on broad-band data used for the continuum flux density can lead to extremely large EW₀ for spurious detections in narrow-band images. This mostly applies to narrow-band selected LAE samples and since the MUSE-Wide and MUSE-Deep samples were constructed using the spectroscopic information of MUSE where we can confirm the presence of an emission line (and thus classify the objects correctly), this danger is lower in our study.

In Table 4 we show an overview of the measured w₀ values for different confidence levels of the classifications (see Sect. 3) and in Table 5 we show the fractions of high EW₀ (EW₀ > 100 Å and EW₀ > 240 Å) for different confidence levels and β parameters. If we exclude objects with a confidence below 2, the characteristic EW₀ values of the EW₀ distribution do not change much and the fraction of high EW₀ (EW₀ > 100 Å and EW₀ > 240 Å) stays almost the same. Using the highest confidence objects reduces our sample to 455 LAEs. The biggest effect is the exclusion of objects with faint (or low S/N) emission lines with potentially smaller EW₀. Although the fraction of HST undetected objects among the highest confidence sample is only ≈12%, compared to around a third for the entire sample, using the high confidence objects results in a similar characteristic EW₀ value w₀ and fraction of high EW₀ objects (both for EW₀ > 100 Å and EW₀ > 240 Å).

The reason why the fraction of HST undetected objects is lower for LAEs with a high confidence is that their Lyman α lines are usually stronger, with a higher flux allowing for a higher classification confidence. A higher Lyman α flux usually comes along with a brighter UV continuum, making the object more likely to be detectable in the HST data as well (see Sect. 6.3 and Fig. 4).

This study demonstrates the existence of LAEs with high EW₀ (above EW₀ > 240 Å) and that they are quite numerous. Notably, one of the highest securely measured EW₀ (with an error that indicates a 3σ confidence) we find is at EW₀ = 588.9 ± 193.4 Å, which is a clear indication of an unusual underlying stellar population (see Sect. 8.3).

6.2. The influence of survey depth on equivalent widths

We now want to understand the influence of cuts in Lyman α luminosity on the measured EW₀ and histograms. Thus, to estimate the effect of different survey depths, we used a statistical experiment, creating a sample of 10 000 objects with UV magnitudes drawn randomly from a luminosity function based on the Schechter parameterisation ϕ^* ln(10) 0.4 × 10^{−0.4(M − M*)(α + 1)}e^{−10−0.4(M − M*)} with parameters according to Bouwens et al. (2015) for objects at redshift z ≈ 3.8 with a characteristic magnitude M^* = −20.88, a normalisation ϕ^* = 1.97 × 10⁻³ Mpc⁻³ and a slope of α = −1.64, in the absolute magnitude range of −23 < M_UV < −17. For the distribution of EW₀ of these objects we assumed w₀ = 80 Å, close to what we measured for our full sample of real LAEs. We drew EW₀ values randomly for the 10 000 objects and obtained the Lyman α line flux by converting the absolute magnitudes from the luminosity function of Bouwens et al. (2015) to continuum flux densities and multiplying that with the EW₀.

The results of this experiment can be seen in Figs. 7 and 8. Lyman α luminosities versus EW₀ are shown in the left panel in Fig. 7. The hard cut on the left in the left panel is caused by the UV continuum magnitude range that we assumed. If we went to fainter magnitudes, it would be possible to populate the area more with higher EW₀ objects for the lower range of Lyman α line fluxes.

Fig. 7. Statistical experiment to understand the cuts of Lyman α luminosity on EW0. Left panel: Lyman α EW0 over Lyman α line luminosity for 10 000 simulated objects (purple dots). The dashed lines indicate different luminosity cuts in the Lyman α line at log10(LLyα[erg s−1]) = [41.5,42,42.5,43]. Right panel: Lyman α EW0 over the UV continuum luminosity. The different colours correspond to the different line luminosity cuts in the left panel.

Fig. 8. Cumulative distributions of EW0 corresponding to the different line flux cuts in Fig. 7. The inset shows the measured w0 values of exponential fits to histograms of EW0 resulting from the luminosity cuts.

Next, we simulated different survey depths by introducing luminosity cuts and re-derived the EW₀ with these new Lyman α luminosities and the corresponding continuum flux densities. The result can be seen in Fig. 8, where we show the different cumulative distributions for the luminosity cuts in Fig. 7. Inset in Fig. 8 are the measured characteristic EW₀ values w₀, which increase with increasing luminosity cut. The same result remains (although with more scatter), if we reduce the number of objects to 1000 and 100, to see if a smaller survey size would change the distribution.

Thus, for a population where the EW₀ distribution is independent of UV luminosity, a survey with shallower depth in Lyman α will preferentially contain a larger fraction of higher EW₀ LAEs than a sample of larger depth. Our result is in agreement with Nilsson et al. (2009), who showed that UV bright selected samples exhibit a lower fraction of high EW₀ LAEs than UV faint selected samples.

6.3. Results: differences between MUSE-Wide and MUSE-Deep

The influence of survey depth on the EW₀ statistics (as described above) can be demonstrated by comparing the measured EW₀ distributions from MUSE-Deep with the MUSE-Wide sample, shown in the histograms in Fig. 9. As can be seen, MUSE-Wide indeed contains a larger fraction of high and extreme EW₀ objects (20% with EW₀ > 240 Å compared to MUSE-Deep having 11% with EW₀ > 240 Å, see Table 3), whereas MUSE-Deep picks up a larger fraction of low-EW₀ objects (62% with EW₀ < 100 Å compared to MUSE-Wide having 47% with EW₀ < 100 Å).

Fig. 9. Histograms of EW0 for MUSE-Wide (green) and MUSE-Deep (blue) in bins of 50 Å normalised to the total number of objects in each survey. The histograms show the distributions including objects with no detectable UV continuum counterpart in HST.

The reason for this discrepancy are the characteristics of the datasets from the two surveys. MUSE-Wide covers a larger survey area, which leads to the discovery of more extreme and thus rare objects with very high EW₀. MUSE-Deep has longer exposure times and thus a lower Lyman α detection flux limit, making it possible to detect fainter Lyman α luminosities. The question is whether the LAEs that are fainter in Lyman α are also fainter in the UV continuum, which would lead to a similar EW₀ distribution between the brighter and the fainter Lyman α objects, keeping in mind that the LAEs were detected based on their Lyman α emission line (although the HST depth is similar in both surveys). It can be seen from the distribution of EW₀ that objects with lower EW₀ are more numerous, especially in the survey going deeper in Lyman α emission.

To investigate the differences between the two surveys, the UV continuum magnitudes (measured here at 1500 Å) and the Lyman α luminosities are shown in Fig. 10, with different colours for MUSE-Wide and -Deep. This plot demonstrates that the combination of a wide and a deep survey is essential to both cover the low-EW₀ and high-EW₀ tails of the LAE EW₀ distribution.

Fig. 10. Logarithmic Lyman α luminosity log10(LLyα[erg s−1]) plotted against the UV continuum magnitude MUV, divided into MUSE-Wide (light green crosses and triangles) and MUSE-Deep (dark green open dots and triangles). The triangles show upper 1σ limits of MUV for HST undetected objects. The diagonal black lines are lines of equal equivalent widths (at the median redshift of z ≈ 4 and using a β value of −1.97) to guide the eye. The grey dashed lines indicate MUV = −18 and log10(LLyα[erg s−1]) = 42. The grey crosses in the top left corner show the median sizes of the errors (for objects with HST counterparts). The black stars show Green Pea galaxies from Izotov et al. (2016a,b, 2018a,b) that were all found to be leaking Lyman continuum.

While objects with strong Lyman α emission are rare (according to Lyman α luminosity functions, see e.g., Drake et al. 2017a; Herenz et al. 2019), objects with faint UV magnitudes are harder to observe, which explains the lack of objects in the upper left part of the figure. However, this is exactly the range in which higher EW₀ objects can be found (lines of equal EW₀ are shown in the plot). Due to the larger survey area of MUSE-Wide, the high Lyman α luminosity part of the plot is mostly populated by MUSE-Wide LAEs, while the range in UV continuum magnitude covered is similar in both samples. MUSE-Deep is detecting more low Lyman α luminosity objects at similar UV continuum flux densities, resulting in many objects with smaller EW₀.

It is interesting to note that some known LyC leaker galaxies at lower redshifts (Green Peas from Izotov et al. 2016a,b, 2018a,b), marked in Fig. 10 as black stars, occupy mostly the MUSE-Wide range of UV continuum magnitude and Lyman α luminosity values. MUSE-Wide LAEs might thus be ideal objects to search for LyC leakage at high redshift.

To test the influence not only of the Lyman α detection limits but also of the UV continuum and survey area, we constructed a cumulative EW₀ distribution in Fig. 11 similar to our statistical experiment in Fig. 8. The fraction of higher EW₀ is higher in MUSE-Wide than in MUSE-Deep, consistent with the discussion above. A Kolmogorov Smirnow (KS) test indicates that the two are likely not drawn from the same underlying distribution (p = 0.0014). Since this plot shows the fraction of objects, the effect of the larger survey area of MUSE-Wide is accounted for. This is because a larger survey area will increase the total number of objects, but not the fraction of objects in certain EW₀ bins. To further test the possible influence of the survey area, we reduced the size of MUSE-Wide to nine fields, to simulate a similar size as MUSE-Deep. This can be seen as the grey lines in Fig. 11, each indicating nine randomly picked MUSE-Wide fields which follow the trend for the full MUSE-Wide survey (albeit with more scatter, since there are fewer objects). Since the distributions of EW₀ for the random MUSE-Wide fields do not match significantly better with the MUSE-Deep distribution (with a median p-value of p = 0.006) than the real MUSE-Wide distribution, the influence of the survey area is likely to be small.

Fig. 11. Cumulative distribution of Lyman α EW0 divided into MUSE-Wide (thick green line) and MUSE-Deep (thick blue line). The lines show the samples including objects with no detected UV continuum in the HST images, each step indicates an individual object. For the objects with UV continuum counterparts, the 1σ error was subtracted from the measured EW0 in order to show the same lower limits as for objects without a continuum counterpart. The grey lines show the cumulative distributions of EW0 for the MUSE-Wide data in an area of a similar size as MUSE-Deep, by picking nine MUSE-Wide fields at random. The black dashed lines indicate EW0 = 100 Å and EW0 = 240 Å.

Next we tried to match the lines for MUSE-Wide and -Deep in the cumulative EW₀ distribution by introducing cuts in the absolute UV magnitude of M_UV < −18 (top panel of Fig. 12) and the Lyman α luminosity log₁₀(L_Lyα[erg s⁻¹]) > 42 (middle and bottom panels). Applying a cut in absolute UV magnitude does not remove but reduces the discrepancy between the two surveys, especially for the high EW₀ where more objects with fainter UV continua would be found. The cut in Lyman α luminosity, shown in the middle and bottom panels of Fig. 12, reduces the discrepancy especially in the low and very high range of EW₀. If we only look at objects with a UV continuum counterpart detected in HST, which effectively also introduces a cut in UV continuum depth, the curves still match well (with the additional Lyman α luminosity cut, see Fig. 12, bottom panel).

Fig. 12. Same plots as Fig. 11 but with a cut in MUV < −18 for the top panel and a cut in Lyman α luminosity of log10(LLyα) > 42 for the middle and lower panel. Both cuts are marked as dashed lines in Fig. 10. In the lower panel, the sample is divided into objects with a UV continuum counterpart (solid lines) and without (dashed lines). As in Fig. 11, the black dashed lines indicate EW0 = 100 Å and EW0 = 240 Å.

This analysis leads to the conclusion that the difference in the EW₀ histogram is mostly due to the deeper Lyman α luminosity detection limit of MUSE-Deep, which leads to more objects with low EW₀. Thus, the construction of the survey has a significant influence on what kind of objects can be observed. Without correcting for these biases, simply looking at the distribution or fraction of high EW₀ (but also just the distribution of EW₀ as a whole) is not enough to understand the real occurrence of such objects. It would be useful to correct the distribution, taking the selection function into account and constructing an EW₀ distribution function, in the style of a luminosity function. We will exemplify such a procedure in a follow-up paper (Kerutt et al., in prep.). However, the wedding-cake approach taken here, using observations with different depth, is already an improvement over previous studies.

7. Connecting equivalent widths to line shape and morphological properties

In this section we study how the Lyman α EW₀ is connected to the shape of the line or of the UV continuum counterpart. In Sect. 8 we go into a more detailed discussion of the properties of the highest EW₀ objects in particular.

7.1. Blue bump fraction and ratio

A common feature of Lyman α radiative transfer models is that they predict the existence of a weaker blue peak in the presence of outflowing gas (e.g., Verhamme et al. 2015). Due to the decreasing IGM transmission for photons to the blue of the rest-frame Lyman α wavelength, the observable fraction of blue bumps at higher redshifts will go down (Laursen et al. 2011 but also Hayes et al. 2021), which is why a redshift cut is useful when discussing fractions of blue bumps. In total, around 33% (324, see Table 6) of our objects below a redshift of z < 4 have a blue bump.

Objects with a blue bump tend to have a slightly larger EW₀ (see Table 6), although with a large spread in values. In fact, the EW₀ values of objects with a blue bump seem to have a larger spread, possibly due to the often smaller S/N of the blue bump, and thus could reach to higher values than the objects without a blue bump. If we assume that most emission lines (more than the approximately one-third we find) have a blue bump which is mostly eaten away by the intervening IGM, then this could explain why we see high EW₀ in cases where the blue bump was not absorbed by the IGM. However, because of the unknown IGM absorption it is not possible in this sample to draw a conclusion about the influence of the presence of a blue bump on a stronger Lyman α emission.

The high fraction of objects with blue bumps matches well with literature results. Kulas et al. (2012) find a similar result of 30% for their LAEs at redshift z ≈ 2 − 3, Trainor et al. (2015) find 40%, while Yamada et al. (2012) even find a fraction as high as 50% in one of the largest overdense regions of LAEs at z = 3.1. In a recent study of lensed LAEs at z ≈ 2 − 3, Cao et al. (2020) also find a blue bump fraction of 30 − 50%. However, most studies omit that the observed fraction of blue bumps in LAE spectra depends on the S/N of the line, with higher detection completeness and blue bump fractions for higher S/N. We can thus assume that most Lyman α spectra have a blue bump that is usually hidden in the noise or absorbed by the IGM. Indeed, Hayes et al. (2021) argue that the decline in the fraction of blue bump to main peak ratio with increasing redshift can be attributed entirely to the rise in the neutral fraction of the IGM. We can see this as well by looking at the fractions of blue bumps in different redshift bins. For the sample of LAEs below a redshift of 4 we find 33% of objects have blue bumps, while the fraction is only 16% (107 LAEs) for objects between 4 < z < 5 and 7% (22 LAEs) for objects above z > 5.

Since the presence and strength of the blue bump are related to radiative transfer processes which can also have an influence on the EW₀, we wanted to test if we see a correlation between the two in our data. In Fig. 13 we compare the ratio between blue bump and total line flux to the rest-frame EW and find no strong correlation between the two, with a p-value of only p = 0.0018 (for a Spearman rank correlation test). It looks like the lowest EW₀ have preferentially large blue bumps with respect to the total line flux, but overall there is a wide spread in fractions. This can also be seen in the histogram, where we divide the sample into high and low EW₀. The p-value for the hypothesis that the two distributions are drawn from the same underlying sample is p = 0.0023 (based on a KS-test), so unlikely but not impossible. Here lower EW₀ can more often be found with stronger blue bumps, which is not what for instance Erb et al. (2014) find for LAEs at z ≈ 2 − 3. They use a slightly different definition, as they do not divide the lines into blue and red peaks, but blue- and redshifted emission along the systemic redshift and find a strong correlation between the fraction the EW₀ of the blue and red parts and EW₀, with a 6.3σ significance. A strong blue bump can indicate the presence of holes or low neutral hydrogen column densities in an outflowing gas cloud. Both could also facilitate the escape of Lyman continuum emission, which would otherwise get absorbed by the neutral hydrogen. Therefore, the strength of the blue bump or the fraction of Lyman α emission that is blueshifted can be used as a proxy for Lyman continuum leakage. It should be noted, though, that for low redshift Lyman continuum leakers, the peak flux ratio does not correlate well with the Lyman continuum escape (see Verhamme et al. 2017), but that the strongest correlation is between Lyman continuum escape and the shift of the blue bump with respect to the systemic redshift (for which we would need additional emission lines to determine it precisely).

Fig. 13. Flux ratio of the peaks shown as a function of the EW0 in Å (with a logarithmic scale). Filled green dots show objects with UV continuum counterparts, empty green dots show objects with lower limits for EW0. The flux ratio is defined as the fraction of the line flux of the blue bump compared to the total line flux (blue bump and red, main peak). The black dots show the median peak ratio (including lower limits) in bins of Δ(log10(EW0)) = 0.2, the width of which is shown as the x-axis error bars. The y-axis error bars show the first and last quartiles. The Spearman rank correlation coefficient and p-value are given in the plot. The grey dashed line marks a ratio of 0.5, meaning above this value, the blue line is stronger than the red peak. The panel on the right shows histograms for the flux ratios of the peaks divided into the half of the objects with higher EW0 (green) and lower EW0 (blue). The p-value based on the Kolmogorov Smirnow (KS) test is shown below the legend of the histogram, indicating the low likelihood that both distributions are taken from a common parent sample.

One caveat of the analysis shown in Fig. 13 is that it is difficult to account for objects without a detectable blue bump, which are not included in this figure. The detectability of a blue bump might itself depend on the ratio between the flux of the blue bump and the total flux, as well as on the S/N and possibly the EW₀, which is why this analysis should be seen as preliminary as these dependencies have to be studied further in a dedicated paper, see Vitte et al. (in prep.).

The median peak ratio (blue bump to total) is f_bb/f_t = 0.237 ± 0.205 (consistent with the average profile found at z ∼ 2, see Matthee et al. 2021), which means the blue bump makes up around 20% of the total flux of the Lyman α line. This only holds for the objects where we could detect a blue bump, but if the ratio is very low and the blue bump thus weak, the line would not be included in this statistic. Therefore the real peak ratio might be lower than what we find here. Notably, there are also 62 of 453 lines (for the full redshift range) where the blue bump is stronger than the red peak, which is a sign of infalling material or an outflow with a low neutral hydrogen column density, which enables the Lyman α photons to escape from the near side of the galaxy (Erb et al. 2010; Alexandroff et al. 2015; Verhamme et al. 2017).

7.2. Connections between line properties

In this subsection we focus on the possible connections between different line properties of the Lyman α line.

7.2.1. Peak separation versus FWHM

As can be seen in Fig. 14, the peak separation and the FWHM show a linear correlation (with a Spearman rank correlation coefficient of 0.5 and a corresponding p-value of 5 × 10⁻³⁰), as lines with a larger FWHM have a higher peak separation. This is because both are tracers of the neutral hydrogen column density in the ISM of the galaxy. We expect this correlation from theoretical radiative transfer models (see Verhamme et al. 2015) and observed by Verhamme et al. (2018) for lower redshift galaxies, as indicated by the white stars. Since a low neutral hydrogen column density lets Lyman α photons escape with less scattering, leading to narrower lines and smaller peak separation, it also allows Lyman continuum emission to escape more easily (e.g., Verhamme et al. 2015; Izotov et al. 2018b). There seems to be a flattening of the relation up to higher peak separations, which might be influenced by small number statistics, though. The peak separation range for possible LyC leakers according to Izotov et al. (2018b) (up to ∼250 km s⁻¹ for a significant Lyman continuum escape fraction above > 10%), marked in grey in Fig. 14, contains 65 objects of our sample (out of 455 objects with a double peak, so ∼12%, see also Vitte et al., in prep.). The reverse conclusion would be, that roughly ∼88% of our objects are not likely to be leaking Lyman continuum emission.

Fig. 14. FWHM of the main (red) peak of the Lyman α line plotted against the peak separation, both in km s−1, colour-coded by the logarithmic EW0. Coloured dots with black circles show objects without UV continuum counterparts. The grey bar in the top left corner shows the median errors. The black dots with error bars are binned median values. The x-axis error bars show the bin widths of 120 km s−1, the y-axis error bars show the first and last quartiles. The grey area is the peak separation range for which Izotov et al. (2018b) show a significant (> 10%) Lyman continuum escape fraction for their sample of low-redshift LyC leakers. The white stars indicate a sample of high-redshift LAEs compiled by Verhamme et al. (2018) (using those objects where the peak separation is given in their Table 1) with measured systemic redshifts. The black line shows their correlation, the dashed line the theoretical line where the peak separation is twice the FWHM. The white dots show a compilation of Green Peas from Yang et al. (2017) (values taken from their Table 3, peak separations are the sum of the individual peak shifts).

The mean peak separation of the full sample of our objects with a visible double peak is 481 ± 244 km s⁻¹ (see Table 7 for an overview of the mean values for the different line properties divided into the high and low EW₀ sample), which can be translated approximately to a shift of the red peak with respect to the systemic redshift of Δv ≈ 241 km s⁻¹ (using the empirical correlation by Verhamme et al. 2018, based on the fact that most Lyman α lines with a double peak are symmetric around the systemic wavelength). In comparison, Shapley et al. (2003) find a shift of Δv ≈ 650 km s⁻¹ for a sample of LBGs, Kulas et al. (2012) find values with a mean of Δv ≈ 370 km s⁻¹ (for half the peak separation), which is closer to our value. Turner et al. (2014) also look at a sample of LBGs and find a difference with respect to systemic redshift of Δv = 220 km s⁻¹ for objects with only Lyman α emission and Δv = 370 km s⁻¹ for those with Lyman α emission and interstellar absorption. Hashimoto et al. (2015) look at LAEs at z ≈ 2.2 and find an even smaller offset between Lyman α and nebular lines of Δv ≈ 174 km s⁻¹, similar to Δv ≈ 180 km s⁻¹ from Song et al. (2014) and Δv = 171  ± 8 km s⁻¹ (difference to systemic redshift) from Muzahid et al. (2020). They both find an anti-correlation between the Lyman α EW₀ and the velocity offset. Since we do not know the systemic redshift of our LAEs, we have to approximate the velocity offset with the peak separation.

7.2.2. Peak separation versus EW₀

Since the peak separation and FWHM are connected to the neutral hydrogen column density, one would expect an anti-correlation between the peak separation and the EW₀ of the lines (see Verhamme et al. 2015), as Lyman α photons scatter in neutral hydrogen, increasing their path length and therefore their probability to be destroyed by dust. This has been shown for a sample of low-z Lyman continuum leaker candidates by Verhamme et al. (2017), but we do not see such a correlation here (see top panel of Fig. 15). In the top right panel in Fig. 15 we show the histogram of peak separation values of the half of the sample having larger EW₀ and the half with smaller EW₀. The half with the higher EW₀ has a narrower distribution, with more objects at low to intermediate peak separations, but using a Kolmogorow-Smirnow-test we can not rule out the null-hypothesis of the same underlying distribution (with a p-value of 0.55).

Fig. 15. Two properties of the Lyman α line plotted against EW0 in Å (in logarithmic scaling). Top panel: peak separation (for objects with a double peaked Lyman α line) in km s−1, bottom panel: FWHM of the main line (for all objects). The symbols are the same as in Fig. 13, the histograms to the right of the plots show again the distribution separated into the half with higher EW0 (in green) and the half with lower EW0 (in blue). The bin width for the histogram for the peak separation is 50 km s−1, the bin width for the histogram for the FWHM is 100 km s−1, the p-value is given based on a KS-test.

7.2.3. FWHM versus EW₀

When it comes to the width of the line, the FWHM, the mean value we find is FWHM = 217.89 ± 102.26 km s⁻¹, which in turn is smaller than what is found by other studies, for instance the average FWHM of Lyman α lines of LBGs at z = 2 − 3 is found to be ≈650 km s⁻¹ by Steidel et al. (2010), and Erb et al. (2010) find an LAE with a FWHM of ≈850 km s⁻¹. The latter study is analysing a single object though and among our sample of LAEs, the widest line has a FWHM that is even higher (901.8 ± 166 km s⁻¹). Another more recent study find FWHM values closer to ours, with FWHM = 212 ± 32 Å for a sample of 3000 LAEs at z ∼ 2 (Hashimoto et al. 2017b).

In the bottom panel of Fig. 15 we find no correlation between the FWHM and the EW₀. However, when we divide our sample again in two parts, one containing the higher EW₀ objects, one the lower, we can see in the histogram in the bottom right panel of Fig. 15 that there is a small difference between the two samples (ruling out the null-hypothesis of the same underlying distribution with a p-value of 5 × 10⁻⁶ using a KS-test). For local LAE analogues, the velocity shift is as small as Δv ≈ 150 km s⁻¹ on average (Verhamme et al. 2017).

7.2.4. Flux ratio versus peak separation

There is an anti-correlation (with a Spearman rank correlation coefficient of −0.41 and a corresponding p-value of 7.9 × 10⁻²⁰) between the contribution of the blue bump to the total line flux and the peak separation in Fig. 16; the peak separation decreases with increasing contribution by the blue bump to the line flux. There are several areas of the plot that have only few objects. One of which is the range of objects with narrow peak separations and small ratios. This could be caused by a bias, as small peak separations make it hard to detect the blue bump due to the line spread function of MUSE and also the intrinsic width of the lines. In addition, small peak ratios are hard to detect as well, which is why this combination might lead to missing objects in the bottom left corner of Fig. 16. There are also no objects with a high peak ratio and a high peak separation.

Fig. 16. Flux ratio of the two peaks plotted as a function of the peak separation, colour-coded by the logarithmic rest-frame EW. As in Fig. 14, the coloured dots with black circles show objects without UV continuum counterparts and the grey bar in the top left corner shows the median errors. The black dots and error bars show the median flux ratio (including lower limits) in peak separation bins of 120 km s−1. The x-axis error bars show the widths of the bins, the y-axis error bars show the first and last quartiles. The grey dashed line shows a ratio of 0.5, meaning above this value the blue bump dominates. The light green stars show Green Pea galaxies form Yang et al. (2016) in a redshift range between z ≈ 0.14 − 0.26.

One reason for this effect could be the IGM absorption of the blue bump, which would decrease its line flux, while the red main peak is shifted out of the resonance frequency. However, the green stars in Fig. 16, showing low redshift Green Peas (Yang et al. 2016), overlap with the measured peak ratios and separations of our sample of high redshift LAEs, highlighting the similarities between these objects. Since the low redshift Green Peas are not affected by absorption in the IGM, the trend can not be explained by this effect. On the one hand, ionised bubbles around LAEs and LBGs have been found even at high redshifts, which could allow a blue bump with a small peak separation to shift out of resonance before reaching the IGM (see Matthee et al. 2018 for a blue bump at redshift z = 6.6 and a peak separation of 220 ± 20 km s⁻¹). The smaller the peak separation, the smaller the ionised bubble needs to be to allow the blue bump to shift enough to survive the IGM. On the other hand, it has been shown in simulations that the larger the peak separation, the more flux of the blue bump will be transmitted through the IGM (see Fig. 2 in Laursen et al. 2011). To learn more about the underlying processes of this anti-correlation, which might be linked to gas kinematics, we need more realistic radiative transfer models.

7.2.5. Asymmetry versus peak separation

The asymmetry of the red part of the line and the peak separation do not show a strong trend (see Fig. 17, with a p-value of 0.04). Only for small values, the peak separation seems to be smaller for higher asymmetries (with a red wing, which is the most common line shape for Lyman α), which matches theoretical models (e.g., Verhamme et al. 2015). On average we can confirm the typically observed line shape of Lyman α in our sample, with an asymmetric line with a red tail and a median asymmetry parameter of a_asym = 0.156.

Fig. 17. Asymmetry of the red peak of the Lyman α line (asymmetry parameter aasym. in Eq. (3)) compared to the peak separation, colour-coded by the logarithmic EW0. As in Fig. 14, coloured dots with black circles show objects without UV continuum counterparts. The grey cross in the top left corner shows the median errors. The black dots and error bars show the median asymmetry (including lower limits) in peak separation bins. The x-axis error bars show the widths of the bins, the y-axis error bars show the first and last quartiles. The asymmetry is defined such that positive values mean the line has a red wing, while negative values mean the line has a blue wing (and zero means no asymmetry, shown as the dashed grey line).

In conclusion, we find that in our sample, LAEs have large peak separations on average, making them similar to values found for LBGs in the literature. Since the distinction between LBGs and LAEs is often an arbitrary cut in EW₀ of ∼ 20 Å which we do not make, this result could be caused by the ability to find faint emission lines in our MUSE data. The widths of our lines are smaller than found in other studies of LAEs, which hints at smaller neutral hydrogen column densities. However, the spread of values is large and we do find extreme values of FWHM ∼ 900 km s⁻¹.

7.3. UV continuum morphology

In this section we explore the connection between Lyman α EW₀ and the UV continuum morphological properties to see if objects with higher EW₀ can be distinguished from the rest of the sample. Since the main driver for the production of Lyman α photons are young, massive stars, an undergoing merger and thus experiencing a boost in star formation, can produce more Lyman α emission and thus higher Lyman α EW₀. In a similar vein, the orientation of the LAE (whether it is seen edge-on or face-on or the Lyman α emission is escaping through a hole in the ISM) could boost the Lyman α photons that are scattered into or out of the line of sight, which in turn influences EW₀ (Verhamme et al. 2012; Behrens et al. 2014; Behrens & Braun 2014; Shibuya et al. 2014a). The orientation results in an elongation of the object, parameterised by the axis ratio which we can measure with Galfit.

Since we fitted the LAEs with a Sérsic profile using Galfit (as explained in Sect. 4.2), we not only obtained the UV magnitudes of the LAEs but also morphological parameters such as the axis ratio, effective radius, position angle, and Sérsic index. While it has been argued based on zoom simulations of high-redshift galaxies that properties such as axis ratio are of limited use when determining the state of the galaxy, such as whether it is merging (Abruzzo et al. 2018), it is interesting to see if the internal properties of the LAE that influence the strength of the Lyman α line, and thus the EW₀, are connected to its morphological properties.

The first morphological property we got from our fits is the number of components of the UV continuum counterpart. As shown in Fig. 2 and discussed in Sect. 4.1, we modelled the counterparts consisting of one (1118 LAEs) or more components (166 LAEs) and found that objects with counterparts of only one component tend to have a larger EW₀ at ${75}_{-44}^{+128}\mathrm{\AA }$ than objects with more components with EW${}_0={49}_{-28}^{+84}\mathrm{\AA }$, however with a considerable spread (errors given as first and last quartiles).

The axis ratio can be used as a proxy for the viewing angle, with rounder objects (higher axis ratio) being viewed more face-on. It has been proposed that face-on galaxies should have a higher Lyman α escape fractions and thus a higher EW₀ (Verhamme et al. 2012; Behrens & Braun 2014), as the Lyman α photons do not have to travel through the dusty disc and can escape more easily. Observationally this has been shown for example by Shibuya et al. (2014b) and Paulino-Afonso et al. (2018), who find a positive correlation between the axis ratio and the rest-frame EW of LAEs in a similar redshift range, however only going up to EW₀ of 300 Å.

Figure 18 shows a comparison between the rest-frame EW₀ and the axis ratio (top panel) and the effective radius R_e (bottom panel). For those comparisons we use the main (brightest) component of the UV continuum counterpart in cases where there were several.

Fig. 18. Two different measures of properties of the UV continuum counterparts of the LAEs (measured in the HST photometry) compared to the rest-frame EW of the Lyman α line. Top panel: axis ratio over EW0, bottom panel: effective radius (in physical kpc) over EW0. As in Fig. 13, the black dots show the median values in bins of Δ(log10(EW0)) = 0.2, indicated by the widths of their error bars in x-direction with the first and last quartiles of the values in each bin indicated by the y-axis error bar. However, for these two plots, objects with a UV continuum counterpart consisting of only one component are shown as filled dots, objects with multiple components are shown in empty dots. For the latter, the brightest component was used to measure the size and axis ratio. The histograms show again the higher EW0 (green) and lower EW0 (blue) half of the sample, the axis ratio in bins of 0.1, and the effective radius in bins of 0.2.

While the axis ratio shows no significant correlation with the rest-frame EW, it is notable (seen from the histogram) that the highest EW₀ objects tend to have either elongated shapes or are very round, while objects with smaller EW₀ values are somewhere in between. In contrast, Shibuya et al. (2014b) find a trend that higher Lyman α rest-frame EW objects have lower UV ellipticity and argue that Lyman α photons can escape more easily from face-on galaxies.

Another noteworthy result is the distribution of the axis ratios (see the histogram in the top right panel of Fig. 18), with a median value of axis ratio = 0.41. The distribution has a maximum and is not flat as seen for local spiral galaxies (e.g., Lambas et al. 1992; Padilla & Strauss 2008).

In Fig. 19 the axis ratio of the LAEs is not correlated with the Lyman α line asymmetry (with a p-value of 0.6). Vivian et al. (2015) have found in their analysis that for more symmetric UV morphological properties, the Lyman α skewness (similar to the asymmetry shown in this work) has a larger scatter, which is not observed in our data. Childs & Stanway (2018) even caution that the measure of the asymmetry or skewness of the line is complicated by sky lines and resolution effects and should rather not be used for interpretations.

Fig. 19. Asymmetry aasym of the Lyman α line (see Sect. 5.2) as a function of the axis ratio of the UV continuum counterpart, similar to Figs. 16 and 17. The black dots show the median values in bins, indicated by the widths of their error bars in x-direction with the first and last quartiles of the values in each bin indicated by the y-axis error bar. The top left cross shows the median errors in both directions, the grey dashed line indicates a symmetric line. Only objects with UV continuum counterparts are shown.

Another morphological property of LAEs is their UV continuum size. In the bottom panel of Fig. 18 there is an anti-correlation (with a p-value of 6.4 × 10⁻¹⁶ for a Spearman rank test) between the effective radius of the UV continuum counterpart and the EW₀. This could be caused by a bias, as faint UV counterparts are more difficult to fit, but will generally have higher EW₀. However, if we include a cut in magnitude of M_UV < −18 to minimise the influence of UV faint objects, the correlation persists. It is also notable that among the objects with the highest EW₀, most are compact with effective radii below ∼1 kpc (physical), comparable for instance with the results in a study of LAE morphological properties by Paulino-Afonso et al. (2018). Looking at the histograms of effective radii divided into the high and low EW₀ half of the sample confirms the anti-correlation, as low EW₀ tend to have larger effective radii. One caveat of this study is that ideally we would need to compare to a sample of galaxies in the same redshift range that do not have Lyman α emission or even show it in absorption.

Looking at the values of the morphological parameters themselves, we find a large variation in the Sérsic index values, indicating a variety of morphological properties in high redshift LAEs, but the median Sérsic index is n = 2 ± 1.8, meaning most LAEs are consistent with an exponential, disc-like profile (as found by e.g., Gronwall et al. 2011 for LAEs at z ≈ 3.1). The median axis ratio is b/a = 0.45 ± 0.21, matching with the value of b/a = 0.45 found by Gronwall et al. (2011), and the distribution is skewed towards more elongated shapes. One should keep in mind, though, that for faint objects, which might be even unresolved in HST, the axis ratio derived from Galfit is influenced by the PSF estimation and the noise (see Gronwall et al. 2011 for a discussion).

8. Discussion of high EW₀ objects

We find a significant number of LAEs with very high EW₀ among our sample, notably 306 LAEs with EW₀ of EW₀ > 240 Å (including both objects with and without counterparts), one of the largest samples of high EW₀ LAEs. In this section we discuss the connections between Lyman α EW₀ and the spectral and morphological properties of the LAEs in order to determine what makes LAEs with high EW₀ (with EW₀ > 240 Å) special.

Results on fractions of LAEs with high EW₀ from previous studies vary widely. Compared to most studies, our result of a fraction of ≈16% of objects with EW₀ > 240 Å is rather high. In a sample of ∼3000 narrow band selected LAEs with redshift z ≈ 2, Hashimoto et al. (2017b) find a sub sample of six spectroscopically confirmed LAEs with an EW₀ of ≈ 200 − 400 Å, which constitutes a lower limit to the fraction of high EW₀ LAEs of ≈0.2%. Similarly, Nilsson et al. (2009) show a fraction of 4% of LAEs at redshift z ∼ 3 with EW₀ > 240 Å and Santos et al. (2020) find only 45 in their SC4K (Sobral et al. 2018) sample of ≈4000 LAEs (≈1.1%) between z ∼ 2 − 6. The SILVERRUSH (Ouchi et al. 2018) survey found a slightly larger fraction, 4% at z ≃ 5.7 and 21% at z ≃ 6.6 (Shibuya et al. 2018). In contrast (at a redshift of z = 4.5), Malhotra & Rhoads (2002) find a fraction of 60% of LAEs exceeding EW₀ > 240 Å in LAEs of the Large Area Lyman Alpha (LALA) survey (Rhoads et al. 2000), however effectively using an IGM correction for the Lyman α line. Using emission line selected LAEs, Adams et al. (2011) find three LAEs (≈3%) at z ∼ 2 − 4 with EW₀ > 240 Å among a sample of 105 galaxies and Shimasaku et al. (2006) find 30 − 40% at z = 5.7. Studies that do not find high EW₀ are for example Gawiser et al. (2006), Gronwall et al. (2007), and Guaita et al. (2011) at redshifts z ≈ 3.1 and z ≈ 2.1. In a recent paper, Maseda et al. (2020) use a sample of UV faint spectroscopically identified LAEs between z = 3 − 4 and stack their photometry to find a median EW₀ = 249 Å with a high ionising photon production efficiency (in combination with Hα measurements). This also shows again that the survey properties are important to take into account when discussing fractions of high EW₀ LAEs, as we argued in Sect. 6.3.

8.1. Connection between high EW₀ and special UV continuum morphological properties

We have seen in the previous Sect. 7.3 that for our sample of LAEs, the half of our sample with larger EW₀ have indistinguishable UV morphological properties to the half with smaller EW₀ when it comes to Sérsic index and axis ratio. As discussed above, the axis ratio could be a proxy for the viewing angle and we could expect that we see higher EW₀ from face-on galaxies. However, we do not find a correlation between axis ratio and EW₀, possibly because the Lyman α photons are indeed escaping primarily perpendicular to the disc’s plane of the LAE, but are scattered back into the line of sight to make them visible from an edge-on view as well. Another explanation is that LAEs do not typically show the same morphological properties as low-redshift spiral galaxies. They are known to be clumpy, often merging systems (even at low redshifts, as can be seen in the LARS sample, e.g., Östlin et al. 2014; Hayes et al. 2014, Messa et al. 2019) and the Lyman α morphological properties do not necessarily follow those of the UV continuum (Östlin et al. 2009; Hayes et al. 2014). Therefore, the interpretation of the axis ratio as a viewing angle might not be accurate for LAEs at these high redshifts, where the galaxies are not well resolved and the axis ratio we measure could be due to irregularities more than the tilt of the LAEs’ galaxy plane.

It has also been proposed that the Lyman α emission strength correlates with the physical size of the UV continuum. Small, compact galaxies could be young, less massive, and less dusty and allow an easier escape of Lyman α emission, such as through their stronger outflows. Our results show that the highest EW₀ objects (above 240 Å) also have smaller effective radii, in agreement with this interpretation. Bond et al. (2012) look at LAEs at redshifts z = 2.1 − 3.1 and find no correlation between the EW₀ and the half-light radius, while for example Marchi et al. (2019) show that UV compact sources have higher Lyman α EW₀ in their sample of star-forming galaxies at z ≈ 3.5. In a sample of LAEs stretching from redshift z ∼ 2 to z ∼ 6, Paulino-Afonso et al. (2018) find a similar trend between the effective radius and the EW₀ as shown in Fig. 18 (bottom panel). Thus we can say that high EW₀ objects (> 240 Å) tend to be compact and small, but small LAEs span a wide range of EW₀.

8.2. Possible causes for high EW₀

From Verhamme et al. (2008) we learn that the FWHM and peak separation of the Lyman α line reflect the neutral hydrogen column density of the ISM in the LAEs. For expanding shell models, line widths of FWHM > 500 km s⁻¹ need neutral hydrogen column densities of N > 10²⁰ cm⁻² and vice versa. Following this model, this suggests that the majority of our objects has small neutral hydrogen column densities below N < 10^19.5 cm⁻², since their FWHM values are on average narrower than ∼230 km s⁻¹ (Verhamme et al. 2008). It has been found in simulations that not only a low neutral hydrogen column density can lead to increased Lyman α EW₀. The production of ionising photons, which are necessary to produce Lyman α emission, is dictated by the internal properties of the galaxy, such as a high star formation rate, a more top-heavy IMF, and a low metallicity (e.g., Schaerer 2003; Raiter et al. 2010). When comparing different IMF models, Raiter et al. (2010) show that for low metallicities (but without invoking PopIII stars yet), EW₀ of up to 500 − 700 Å can be reached, even up to 1500 Å when assuming top-heavy IMF models and a departure from Case B. The number of ionising photons produced in such simulations also depends on the stellar population synthesis models that were used. Models such as BPASS (Eldridge & Stanway 2009, which however assumes a binary fraction of 100%) can produce much higher numbers of ionising photons and thus also higher Lyman α EW₀, up to intrinsic EW₀ of over 1000 Å (see Garel et al. 2021). The hypothesis that indeed such unusual stellar populations are responsible for the high EW₀ could be investigated with further observations, for example probing young stellar ages with Balmer lines, massive stars with N V P-Cygni profiles or measuring the metallicity (e.g., using R23 = (O II + O III)/Hβ), which is expected to be low for objects with high EW₀ (EW₀ > 240 Å).

Directional boosting of Lyman α emission due to a clumpy, multi-phase medium as proposed by Neufeld (1991) could also explain high EW₀ (EW₀ > 240 Å). It has been found to work only under special conditions, though, if the ISM is very clumpy, with clumps over five times more dense than typical, very dusty, and almost static with outflow velocities less than 10 km s⁻¹ (Laursen et al. 2013; Gronke & Dijkstra 2014; Duval et al. 2014; Blaizot et al., in prep.).

To test this prediction, we can imagine that a static ISM would result in double peaked lines where both peaks contribute equally to the total flux, which means we could expect a trend that for a flux ratio of bb/total flux ≈50% the equivalent widths could be boosted, which we can not confirm (see Fig. 13). However, we do not take into account the possible absorption of the blue part of the line by the IGM, which could be a reason why we do not see a correlation here. However, we can see in Fig. 16 that local galaxies which are not affected by IGM absorption seem to follow a similar trend as our high-redshift LAEs, when it comes to the flux ratio, which we infer from the comparison between the flux ratio and the peak separation which fall on the same correlation for our high-redshift LAEs and the local Green Peas from Yang et al. (2016). If the IGM absorption of the blue bump was significant, the high-redshift LAEs would fall below the correlation of the local Green Peas, as the flux ratio would be reduced. Therefore we can believe that indeed we do not see a boosted EW₀ from a static ISM in Fig. 13.

Instead of explaining high EW₀ with a multi-phase, possibly static ISM, Laursen et al. (2013) suggest alternative production mechanisms for large EW₀ values of Lyman α photons, such as collisional emission (e.g., Dijkstra 2017; Faucher-Giguère et al. 2017) and fluorescence (e.g., Cantalupo et al. 2005; Kollmeier et al. 2010; Dijkstra 2017). However, directional boosting of Lyman α EW₀ can work (although with a large directional variability) when it is related to outflows caused by starbursts in combination with higher UV continuum absorption (Smith et al. 2019), resulting in narrower, single peaked lines.

A small peak separation, possibly indicative of a lower neutral hydrogen column density (according to radiative transfer simulations), is related to the escape of Lyman continuum emission, which can produce an ionised region, facilitating the escape of blue bump photons. There could also be a selection bias in this plot, since weak blue bumps, resulting in a small peak ratio, that are also close to the main red peak (with a small peak separation), might not be easily detected due to instrumental limitations (see Matthee et al. 2018 for a discussion).

In that sense LAEs with high EW₀ (EW₀ > 240 Å) are also interesting for the study and search for Lyman continuum leakers, which are responsible for providing the ionising emission at the epoch of reionisation. LAEs with high EW₀ have been found to have elevated ξ_ion (ionising photon production efficiency) values, lower metallicities, and younger ages (Maseda et al. 2020). We have seen for example in Fig. 10 that the Green Pea galaxies from Izotov et al. (2016a,b, 2018a,b), all of which are leaking Lyman continuum, fall in a region of higher Lyman α EW₀. Since a higher Lyman α EW₀ not only indicates a strong production of Lyman continuum photons, but also the possibility of a high escape fraction of Lyman α emission (especially considering that the half of our sample with higher and lower EW₀ have potentially different line widths, see Fig. 15 and the discussion in Sect. 7.2), it has been proposed to use high Lyman α EW₀ as an indirect tracer of the Lyman continuum (see Verhamme et al. 2015, 2017; Marchi et al. 2018; Steidel et al. 2018; Fletcher et al. 2019, but see also Bian & Fan 2020 who contest this). This makes our sample of 307 objects with EW₀ > 240 Å all the more interesting for possible follow-up observations to confirm the presence of Lyman continuum leakage.

8.3. The highest EW₀ LAE

Among our sample, the highest robust EW₀ is EW₀ = 588.9 ± 193.4 Å (which is the highest EW₀ that is measured with a 3σ certainty), with the ID 123007091. The highest lower limit EW₀, where no HST counterpart could be identified, is EW₀ = 4464 Å. Again, the choice of β influences the EW₀. For β = −1.57 (which is the upper quartile of the distribution of measured β values, see Sect. 6.1) the object with the highest robust measurement would have EW₀ = 689.26 ± 260.5 Å, while for β = −2.29 (the lowest quartile) it would have EW₀ = 519.3 ± 153.0 Å. The Lyman α luminosity is log₁₀(L_Lyα[erg s⁻¹]) = 42.56 ± 0.26 and with a UV magnitude of M_UV = −17.2 ± 0.3 it is rather faint.

Kashikawa et al. (2012) find an even higher EW₀ of ${\mathrm{EW}}_0={872}_{-298}^{+844}\mathrm{\AA }$ at a redshift of z = 6.5, which could hint at a candidate population III galaxy. Their measured a Lyman α luminosity of (1.07 ± 0.30)×10⁴³ erg s⁻¹. Assuming a fixed β = 2, they get a UV continuum flux density at the Lyman α wavelength of f_cont = (7.12 ± 3.41) × 10⁻²¹ erg s⁻¹ cm⁻²/Å, compared to f_cont = (15.22 ± 4.68) × 10⁻²¹ erg s⁻¹ cm⁻²/Å for our object.

However, they apply an IGM correction and assume that around half of the Lyman α emission line flux was absorbed by the IGM. Their uncorrected value is at ${\mathrm{EW}}_0={436}_{-149}^{+422}\mathrm{\AA }$, lower than our highest value and with a large error, indicating that it could potentially be much higher.

It is interesting to look more closely at the properties of our object 123007091, which has a redshift of z = 3.26 (see Fig. 20 for the Lyman α halo measurements and Fig. 21 for the UV continuum). The spectrum shows no significant UV emission lines apart from Lyman α. The Lyman α emission line has a tentative close blue bump at a peak separation of 197 km s⁻¹ and a FWHM of the main, red peak of FWHM = 178 ± 110 km s⁻¹ (corrected for the MUSE LSF). Both the line separation and the line width are smaller than the average of the full sample (shown in Table 7), although the FWHM is within the one sigma range of the average. Still, based on the peak separation, this indicates that the Lyman α emission did not scatter much to escape the galaxy, possibly due to a low neutral hydrogen column density or clear paths in the ISM in our direction of observation, making it relatively easy for the Lyman α emission and possibly also the Lyman continuum emission to escape. This object has a spatially compact Lyman α emission with no measurable extended halo beyond the UV continuum (see Fig. 20), which supports the assumption that Lyman α (and also Lyman continuum) escapes more easily from objects with compact Lyman α emission, producing a higher EW₀.

Fig. 20. Measurements of the Lyman α halo of object 123007091 as in Leclercq et al. (2017). The left panel shows the MUSE Lyman α narrow-band image, smoothed slightly with the PSF FWHM Gaussian kernel, indicated by the white hatched circle in the lower right corner. The black lines show the 2σ surface brightness limit, the white dashed circle shows the extraction aperture for the spectrum of the Lyman α emission line. The latter is shown in the middle panel, with vertical lines indicating the width of the narrow-band image. The right panel shows the extent of the Lyman α halo in blue and the continuum component in orange, as well as the PSF (red dotted line).

Fig. 21. HST cutouts of the bands ACS F606W, F775W, F814W and WFC3 F125W and F160W. For the object shown in Fig. 20. As in Fig. 1, the cross is centred on the UV continuum position, the circle is centred on the Lyman α position from LSDCat, and has have a diameter of $0\stackrel{″}{.}5$.

9. Summary and conclusion

In this paper we presented the derived EW₀ for a sample of ∼2000 LAEs with homogeneous measurements of their properties. The LAE sample was constructed using a consistent procedure for the detection and fitting of emission lines (without applying an IGM correction) in data from the MUSE-Wide and MUSE-Deep GTO surveys. We showed EW₀ histograms (Figs. 6 and 9) and an analysis of the Lyman α and morphological properties of the LAEs. Furthermore we connected the measured EW₀ to other spectral and UV continuum morphological properties and discussed the discovery of numerous high EW₀ objects (with EW₀ > 240 Å). The main results of our work are the following:

• A third of our LAEs do not have UV continuum counterparts visible in the HST broad-band data. Of the remaining objects, 13% have multi-component counterparts, hinting at merging, clumpy morphology or star formation (see Figs. 1 and 2).

• A fraction of ≈16% (306 objects) of the full sample of our LAEs have EW₀ > 240 Å, which is comparable to other studies in the literature, especially when considering Lyman α EW₀ studies that, like us, do not apply an IGM correction (see Sect. 8). Due to our careful measurements of the EW₀ values, we establish firmly the existence of a population of galaxies with very high EW₀ (with EW₀ > 240 Å and even much above, see e.g., Figs. 11, 13 and 15).

• In addition we show that such fractions of high EW₀ > 240 Å are influenced by the Lyman α line flux depth of the survey, since we find ≈20% for the MUSE-Wide LAEs and ≈11% for the ones from MUSE-Deep.

• There is a discrepancy between the EW₀ distributions of the MUSE-Wide and MUSE-Deep survey LAEs (Fig. 9), with the deeper data having a smaller fraction of high EW₀ > 240 Å objects than the shallower survey with a larger survey area. We propose as an explanation the ability to detect fainter Lyman α emission in the deeper data. Thus a deeper Lyman α survey at fixed UV continuum depth unveils a larger fraction of the galaxy population which typically has a lower EW₀ (e.g., Shapley et al. 2003 find only an EW₀ of 14.3 Å in their stack of LBGs at z ∼ 3).

• Smaller Lyman α luminosities allow to detect the numerous UV-bright objects with small EW₀ rather than the rare UV faint extreme EW₀ objects.

• We find an object with an extremely high EW₀ of EW₀ = 588.93 ± 193.4 Å, hinting at unusual stellar populations. Among the objects with high EW₀, there is a large spread in emission line shapes (as exemplified in Fig. 5), which implies a large diversity of the influence of radiative transfer processes and the interplay between gas dynamics and dust content. Therefore, the underlying cause of high EW₀ (EW₀ > 240 Å) is probably rooted less in extreme transfer processes, but more in the properties of the stellar populations with respect to age, IMF, and metallicity.

• We confirm the correlation reported in the literature between the effective radius of the UV continuum counterpart and the Lyman α EW₀ (bottom panel of Fig. 18), but do not find a correlation between Lyman α EW₀ and axis ratio (top panel of Fig. 18), which should be a proxy for the viewing angle. We conclude that compact galaxies make it easier for Lyman α to escape, while the axis ratio of the UV continuum counterpart might not be correlated to the viewing angle, but influenced by irregular shapes or star formation regions.

• Objects belonging to the half of the sample with higher EW₀ have potentially different line width properties than the half with lower EW₀ (Fig. 15).

• In general we find that objects with large EW₀ (if we consider the half of our sample with higher EW₀) can have a wide range of spectral line shapes and UV morphological properties (but preferentially small sizes), which means underlying properties of the LAEs can also be varied and there is no only one driving force for producing high EW₀ values.

As we have discussed in this paper, the EW₀ distribution and the fraction of high EW₀ (> 240 Å) objects depend heavily on the survey design that the sample of LAEs was taken from. The observation depth and the area of the survey determine the number of Lyman α faint and bright sources that can be detected. A wedding-cake strategy as used here with the MUSE-Wide and -Deep surveys allows us to find a wider range of LAE properties. A deep survey will unveil more faint LAEs with lower EW₀, which is useful for constraining the faint-end of the LF (e.g., Drake et al. 2017a), while a larger survey area at lower exposure times will uncover more rare objects with large EW₀ (given that the UV continuum depth stays the same). Therefore, to understand the true nature of the LAE population with respect to EW₀, one has to take into account the survey properties. For this it would be ideal to construct an EW₀ distribution function giving the number of LAEs in bins of EW₀ per volume, similar to a luminosity function, thus correcting for the selection function in the two surveys. We will explore this possibility in an upcoming paper.

It is also interesting to study the reasons for the high EW₀ of EW₀ > 240 Å which we find in 307 LAEs. They are ideal candidates to observe with the upcoming JWST instrument to better understand their properties, for example using Hα and Hβ to infer the produced Lyman α luminosity and thus to get the Lyman α escape fraction and learn about the nature of the stellar population. To explore the line shape properties in full it will also be useful to obtain systemic redshifts from JWST. We would also be able to search for rest-frame nebular lines and obtain metallicities, stellar ages, and dust content to see if we can explain their astonishingly high Lyman α EW₀ by extreme ages and metallicities or if we need even more exotic scenarios.

Acknowledgments

The authors wish to thank the anonymous referee for the thorough comments and insight which have helped to significantly improve the quality of the manuscript. We would further like to thank Peter Weilbacher for carefully reading this manuscript and providing suggestions to greatly improve the legibility and conciseness. We also thank Pascal Oesch, Daniel Schaerer, and Moupiya Maji for helpful discussions. Special thanks go to T. Dumont for a fruitful discussion about cinnamon stars. No (LSD)cats were harmed in the making of this paper. The plots in this paper were created using Matplotlib (Hunter 2007) and we also used the Python packages Numpy (van der Walt et al. 2011) and Astropy (Astropy Collaboration 2013). J.K., A.V., H.K., and F.L. are supported by the SNF Professorship PP00P2_176808. A.V., V.M., and T.G. are supported by the ERC Starting Grant 757258 ‘TRIPLE’. J.K., E.C.H., T.U., K.S., and L.W. acknowledge funding from the Competitive Fund of the Leibniz Association through grants SAW-2013-AIP-4 and SAW-2015-AIP-2. H.K. acknowledges support from the Swiss Government Excellence Scholarships and Japan Society for the Promotion of Science (JSPS) Overseas Research Fellowship. T.H. was supported by Leading Initiative for Excellent Young Researchers, MEXT, Japan (HJH02007). We also want to thank the staff at the VLT during the observations of the GTO data.

References

Appendix A: Assessment of the influence of halo sizes

It has been shown in many studies that high redshift Lyman α emitters usually have extended Lyman α halos (Steidel et al. 2011; Momose et al. 2014, 2016; Wisotzki et al. 2016; Leclercq et al. 2017; Wisotzki et al. 2018; Leclercq et al. 2020) that can carry the majority of the emitted Lyman α flux (≈65% as shown by Leclercq et al. 2017). Since the EW₀ depends on correct measurements of the full Lyman α flux from a galaxy, we need to understand if we are capturing the contribution of the halo emission reasonably well. Although there exist measurements of the line flux emission that include the Lyman α halo for some of the objects in this paper, not all are bright enough to perform a curve-of-growth measurement as was done for example by Leclercq et al. (2017), Drake et al. (2017b,a), and Hashimoto et al. (2017a). In order to use a consistent method to obtain the emission line fluxes for all objects in the sample of this paper (both the objects in the MUSE-Wide fields as well as the MUSE-Deep fields), we compared the line flux measurements of LSDCat (which are available for all objects in our sample) to the selected sample of objects from Leclercq et al. (2017), where the halo and UV continuum contributions of the Lyman α halo were analysed separately to obtain the full line flux as well as the halo sizes and the contribution of the halo to the line flux. In Fig. A.1 we show a comparison between the line flux values from Hashimoto et al. (2017a) and Leclercq et al. (2017) that specifically include the halo contribution and line flux measurements within three Kron radii from LSDCat (the same method that was used for the MUSE-Wide fields by Urrutia et al. 2019). It can be seen that the measurements from LSDCat match well for most objects, for both the Hashimoto et al. (2017a) and Leclercq et al. (2017) line flux values. However, there is a slight (and expected) trend that for objects with a large halo scale length (in green and yellow in Fig. A.1), the line flux measurements from LSDCat are smaller and thus do not capture the full halo contribution. Nevertheless, except for those few extreme objects, the line flux measurements match well, which is why we use the LSDCat measurements throughout this paper in order to keep the method consistent. One should keep in mind, though, that the EW₀ estimates are on average slightly underestimated for objects with large halos.

Fig. A.1. Comparison of line flux measurements from different catalogues and methods. The Lyman α line flux from Hashimoto et al. (2017a) (top panel, on the x-axis) and the Lyman α line flux from Leclercq et al. (2017) (bottom panel) are compared to the measurements from LSDCat (y-axis, in the same fields using the same observational depths). The one-to-one relation is shown as the dashed line and the colour-coding and sizes of the dots correspond to the halo scale lengths rshalo from Leclercq et al. (2017). The line flux measurement from Hashimoto et al. (2017a) is based on a curve-of-growth measurement from narrow-band images as described by Drake et al. (2017a,b), similar to Leclercq et al. (2017) who use an additional two-dimensional decomposition for the halo and continuum components of the LAEs.

Appendix B: Table of measurements

A full table containing the measurements shown in this study can be found in the online version of this paper. An overview of the different columns and an example of 34 objects is shown in Table B.1.

All Tables

All Figures

Fig. 1. Examples of LAEs with different kinds of UV continuum counterparts, shown in 2″ × 2″ cutouts of the filter band ACS F814W. The white crosses are centred on the UV continuum position, the white circles are centred on the Lyman α position from LSDCat and have a diameter of $0\stackrel{″}{.}5$. Left image: object with multiple components or clumps in the counterpart. Middle image: object without a significant detection in the filter band ACS F814W (the centre of the cross is set on a noise peak/bright pixel). Right image: counterpart with only one component, but slightly offset with respect to the Lyman α position.

In the text

Fig. 2. Redshift histogram of the full sample of LAEs, split into groups according to their UV continuum counterparts in HST, in bins of Δz = 0.25. The dark blue histogram shows objects with UV counterparts which are not resolved into multiple components. The turquoise histogram shows objects without a visible UV counterpart in the available HST imaging and the light green histogram shows objects with UV counterparts consisting of more than one component. In these cases, the HST counterpart could not be fitted with a single Sérsic profile, as it was clumpy or asymmetric.

In the text

Fig. 3. Screenshot of the Galfit wrapper used for fitting and measuring the UV continuum in different bands. The left image shows a cutout of the data in the HST filter band ACS F814W, the middle image shows the (in this case) two models for the object and the right image shows the residuals. Clicking in the right panel will allow the user to add a Galfit model for an additional component, which creates a green circle to mark the place. Below the panels are morphological measurements as well as magnitude and position, each of which can be fixed if needed. In the example, the object was fitted with two Sérsic models, each of which can be selected by clicking on the small green circles. The check boxes on the right allow the user to asses whether the object is clumpy and if there are multiple components in the counterpart or no counterpart at all, as well whether to include the model in the object.

In the text

	Fig. 4. Histogram showing the logarithmic Lyman α luminosity of the LAEs in the MUSE-Wide (light green, labelled MW) and MUSE-Deep surveys (dark green, labelled MD). The samples are divided into objects with UV continuum counterparts and without counterparts. The number of objects per sample are stacked, which means every histogram bar gives the total number of objects in that luminosity bin.
In the text

Fig. 5. Overview of the Lyman α emission lines of the 35 LAEs with the highest measured EW0, sorted by EW0 in descending order (left to right, top to bottom). We note that the y-axis shows the flux normalised to the maximum of each peak for better comparison. The vertical dashed lines indicate the positions of the Lyman α lines and the blue bump if there is one. The grey shaded area shows the standard deviation from the MUSE spectra. The shown objects are taken both from the MUSE-Wide and -Deep surveys and have been continuum subtracted.

In the text

Fig. 6. Histograms of the Lyman α EW0 of the full sample, including both MUSE-Wide and MUSE-Deep objects. The blue histogram shows the objects with a secure measurement (excluding lower limits), the green bars on top show the objects with a limiting EW0. The exponential fits N = N0 exp(−EW0/w0) are shown as the dashed lines, blue for the objects with secure measurements and light green for the full sample including lower limits. The bin width in EW0 is 50 Å. For the measurement of the characteristic EW0w0 (given below the legend and in Table 4) the smallest EW0 bin was excluded as it is likely incomplete (as Lyman α lines with small or zero line fluxes are hard to detect and preferentially result in small EW0). The shaded areas around the dashed lines indicate the exponential fits for the EW0 distributions using UV continuum slopes of β = −1.57 and −2.29, the first and last quartiles of the distribution of measured β values.

In the text

Fig. 7. Statistical experiment to understand the cuts of Lyman α luminosity on EW0. Left panel: Lyman α EW0 over Lyman α line luminosity for 10 000 simulated objects (purple dots). The dashed lines indicate different luminosity cuts in the Lyman α line at log10(LLyα[erg s−1]) = [41.5,42,42.5,43]. Right panel: Lyman α EW0 over the UV continuum luminosity. The different colours correspond to the different line luminosity cuts in the left panel.

In the text

	Fig. 8. Cumulative distributions of EW0 corresponding to the different line flux cuts in Fig. 7. The inset shows the measured w0 values of exponential fits to histograms of EW0 resulting from the luminosity cuts.
In the text

	Fig. 9. Histograms of EW0 for MUSE-Wide (green) and MUSE-Deep (blue) in bins of 50 Å normalised to the total number of objects in each survey. The histograms show the distributions including objects with no detectable UV continuum counterpart in HST.
In the text

Fig. 10. Logarithmic Lyman α luminosity log10(LLyα[erg s−1]) plotted against the UV continuum magnitude MUV, divided into MUSE-Wide (light green crosses and triangles) and MUSE-Deep (dark green open dots and triangles). The triangles show upper 1σ limits of MUV for HST undetected objects. The diagonal black lines are lines of equal equivalent widths (at the median redshift of z ≈ 4 and using a β value of −1.97) to guide the eye. The grey dashed lines indicate MUV = −18 and log10(LLyα[erg s−1]) = 42. The grey crosses in the top left corner show the median sizes of the errors (for objects with HST counterparts). The black stars show Green Pea galaxies from Izotov et al. (2016a,b, 2018a,b) that were all found to be leaking Lyman continuum.

In the text

Fig. 11. Cumulative distribution of Lyman α EW0 divided into MUSE-Wide (thick green line) and MUSE-Deep (thick blue line). The lines show the samples including objects with no detected UV continuum in the HST images, each step indicates an individual object. For the objects with UV continuum counterparts, the 1σ error was subtracted from the measured EW0 in order to show the same lower limits as for objects without a continuum counterpart. The grey lines show the cumulative distributions of EW0 for the MUSE-Wide data in an area of a similar size as MUSE-Deep, by picking nine MUSE-Wide fields at random. The black dashed lines indicate EW0 = 100 Å and EW0 = 240 Å.

In the text

Fig. 12. Same plots as Fig. 11 but with a cut in MUV < −18 for the top panel and a cut in Lyman α luminosity of log10(LLyα) > 42 for the middle and lower panel. Both cuts are marked as dashed lines in Fig. 10. In the lower panel, the sample is divided into objects with a UV continuum counterpart (solid lines) and without (dashed lines). As in Fig. 11, the black dashed lines indicate EW0 = 100 Å and EW0 = 240 Å.

In the text

Fig. 13. Flux ratio of the peaks shown as a function of the EW0 in Å (with a logarithmic scale). Filled green dots show objects with UV continuum counterparts, empty green dots show objects with lower limits for EW0. The flux ratio is defined as the fraction of the line flux of the blue bump compared to the total line flux (blue bump and red, main peak). The black dots show the median peak ratio (including lower limits) in bins of Δ(log10(EW0)) = 0.2, the width of which is shown as the x-axis error bars. The y-axis error bars show the first and last quartiles. The Spearman rank correlation coefficient and p-value are given in the plot. The grey dashed line marks a ratio of 0.5, meaning above this value, the blue line is stronger than the red peak. The panel on the right shows histograms for the flux ratios of the peaks divided into the half of the objects with higher EW0 (green) and lower EW0 (blue). The p-value based on the Kolmogorov Smirnow (KS) test is shown below the legend of the histogram, indicating the low likelihood that both distributions are taken from a common parent sample.

In the text

Fig. 14. FWHM of the main (red) peak of the Lyman α line plotted against the peak separation, both in km s−1, colour-coded by the logarithmic EW0. Coloured dots with black circles show objects without UV continuum counterparts. The grey bar in the top left corner shows the median errors. The black dots with error bars are binned median values. The x-axis error bars show the bin widths of 120 km s−1, the y-axis error bars show the first and last quartiles. The grey area is the peak separation range for which Izotov et al. (2018b) show a significant (> 10%) Lyman continuum escape fraction for their sample of low-redshift LyC leakers. The white stars indicate a sample of high-redshift LAEs compiled by Verhamme et al. (2018) (using those objects where the peak separation is given in their Table 1) with measured systemic redshifts. The black line shows their correlation, the dashed line the theoretical line where the peak separation is twice the FWHM. The white dots show a compilation of Green Peas from Yang et al. (2017) (values taken from their Table 3, peak separations are the sum of the individual peak shifts).

In the text

Fig. 15. Two properties of the Lyman α line plotted against EW0 in Å (in logarithmic scaling). Top panel: peak separation (for objects with a double peaked Lyman α line) in km s−1, bottom panel: FWHM of the main line (for all objects). The symbols are the same as in Fig. 13, the histograms to the right of the plots show again the distribution separated into the half with higher EW0 (in green) and the half with lower EW0 (in blue). The bin width for the histogram for the peak separation is 50 km s−1, the bin width for the histogram for the FWHM is 100 km s−1, the p-value is given based on a KS-test.

In the text

Fig. 16. Flux ratio of the two peaks plotted as a function of the peak separation, colour-coded by the logarithmic rest-frame EW. As in Fig. 14, the coloured dots with black circles show objects without UV continuum counterparts and the grey bar in the top left corner shows the median errors. The black dots and error bars show the median flux ratio (including lower limits) in peak separation bins of 120 km s−1. The x-axis error bars show the widths of the bins, the y-axis error bars show the first and last quartiles. The grey dashed line shows a ratio of 0.5, meaning above this value the blue bump dominates. The light green stars show Green Pea galaxies form Yang et al. (2016) in a redshift range between z ≈ 0.14 − 0.26.

In the text

Fig. 17. Asymmetry of the red peak of the Lyman α line (asymmetry parameter aasym. in Eq. (3)) compared to the peak separation, colour-coded by the logarithmic EW0. As in Fig. 14, coloured dots with black circles show objects without UV continuum counterparts. The grey cross in the top left corner shows the median errors. The black dots and error bars show the median asymmetry (including lower limits) in peak separation bins. The x-axis error bars show the widths of the bins, the y-axis error bars show the first and last quartiles. The asymmetry is defined such that positive values mean the line has a red wing, while negative values mean the line has a blue wing (and zero means no asymmetry, shown as the dashed grey line).

In the text

Fig. 18. Two different measures of properties of the UV continuum counterparts of the LAEs (measured in the HST photometry) compared to the rest-frame EW of the Lyman α line. Top panel: axis ratio over EW0, bottom panel: effective radius (in physical kpc) over EW0. As in Fig. 13, the black dots show the median values in bins of Δ(log10(EW0)) = 0.2, indicated by the widths of their error bars in x-direction with the first and last quartiles of the values in each bin indicated by the y-axis error bar. However, for these two plots, objects with a UV continuum counterpart consisting of only one component are shown as filled dots, objects with multiple components are shown in empty dots. For the latter, the brightest component was used to measure the size and axis ratio. The histograms show again the higher EW0 (green) and lower EW0 (blue) half of the sample, the axis ratio in bins of 0.1, and the effective radius in bins of 0.2.

In the text

Fig. 19. Asymmetry aasym of the Lyman α line (see Sect. 5.2) as a function of the axis ratio of the UV continuum counterpart, similar to Figs. 16 and 17. The black dots show the median values in bins, indicated by the widths of their error bars in x-direction with the first and last quartiles of the values in each bin indicated by the y-axis error bar. The top left cross shows the median errors in both directions, the grey dashed line indicates a symmetric line. Only objects with UV continuum counterparts are shown.

In the text

Fig. 20. Measurements of the Lyman α halo of object 123007091 as in Leclercq et al. (2017). The left panel shows the MUSE Lyman α narrow-band image, smoothed slightly with the PSF FWHM Gaussian kernel, indicated by the white hatched circle in the lower right corner. The black lines show the 2σ surface brightness limit, the white dashed circle shows the extraction aperture for the spectrum of the Lyman α emission line. The latter is shown in the middle panel, with vertical lines indicating the width of the narrow-band image. The right panel shows the extent of the Lyman α halo in blue and the continuum component in orange, as well as the PSF (red dotted line).

In the text

	Fig. 21. HST cutouts of the bands ACS F606W, F775W, F814W and WFC3 F125W and F160W. For the object shown in Fig. 20. As in Fig. 1, the cross is centred on the UV continuum position, the circle is centred on the Lyman α position from LSDCat, and has have a diameter of $0\stackrel{″}{.}5$.
In the text

Fig. A.1. Comparison of line flux measurements from different catalogues and methods. The Lyman α line flux from Hashimoto et al. (2017a) (top panel, on the x-axis) and the Lyman α line flux from Leclercq et al. (2017) (bottom panel) are compared to the measurements from LSDCat (y-axis, in the same fields using the same observational depths). The one-to-one relation is shown as the dashed line and the colour-coding and sizes of the dots correspond to the halo scale lengths rshalo from Leclercq et al. (2017). The line flux measurement from Hashimoto et al. (2017a) is based on a curve-of-growth measurement from narrow-band images as described by Drake et al. (2017a,b), similar to Leclercq et al. (2017) who use an additional two-dimensional decomposition for the halo and continuum components of the LAEs.

In the text