The MUSE Hubble Ultra Deep Field Survey XIV. The evolution of the Lya emitter fraction from z=3 to z=6
Haruka Kusakabe Jeremy Blaizot Thibault Garel Anne Verhamme Roland Bacon Johan Richard Takuya Hashimoto Hanae Inami Simon Conseil Bruno Guiderdoni Alyssa B. Drake Edmund Christian Herenz Joop Schaye Pascal Oesch Jorryt Matthee Raffaella Anna Marino Kasper Borello Schmidt Roser Pello Michael Maseda Floriane Leclercq Josephine Kerutt Guillaume Mahler
astro-ph.GA astro-ph.CO
 


   
    XIV. The evolution of the Lyα emitter fraction from z=3 to z=6

   

   

Observatoire de Genève, Université de Genève, 51 chemin de Pégase, 1290 Versoix, Switzerland

Univ. Lyon, Univ. Lyon1, Ens de Lyon, CNRS, Centre de Recherche Astrophysique de Lyon UMR5574, F-69230, Saint-Genis-Laval, France


Tomonaga Center for the History of the Universe (TCHoU), Faculty of Pure and Applied Sciences, University of Tsukuba,  Tsukuba, Ibaraki 305-8571, Japan  

Hiroshima Astrophysical Science Center, Hiroshima University 1-3-1 Kagamiyama, Higashi-Hiroshima, Hiroshima 739-8526, Japan,
Gemini Observatory/NSF’s OIR Lab, Casilla 603, La Serena, Chile,

Max-Planck-Institut für Astronomie, Königstuhl 17, D-69117 Heidelberg, Germany

ESO Vitacura, Alonso de Córdova 3107,Vitacura, Casilla 19001, Santiago de Chile, Chile 

Leiden Observatory, Leiden University, PO Box 9513, 2300 RA Leiden, The Netherlands 

Department of Physics, ETH Zürich, Wolfgang-Pauli-Strasse 27, 8093 
Zürich, Switzerland

Leibniz-Institut für Astrophysik Potsdam, An der Sternwarte 16 
14482 Potsdam, Germany

Aix Marseille Université, CNRS, CNES, LAM (Laboratoire d’Astrophysique de Marseille), UMR 7326, 13388, Marseille, France

Department of Astronomy, University of Michigan, 1085 South University Ave, Ann Arbor, MI 48109, USA 










 
  
  
   The Lyα emitter (LAE) fraction, X_ LAE, is a potentially powerful probe of the evolution of the intergalactic neutral hydrogen gas fraction. However, uncertainties in the measurement of X_ LAE are still debated.
  
  
   Thanks to deep data obtained with the integral field spectrograph MUSE (Multi-Unit Spectroscopic Explorer), we can measure the evolution of the LAE fraction homogeneously over a wide redshift range of z≈3–6 for UV-faint galaxies (down to UV magnitudes of M_1500≈-17.75). This is significantly fainter than in former studies (M_1500≤-18.75), and allows us to probe the bulk of the population of high-redshift star-forming galaxies. 
  
   We construct a UV-complete photometric-redshift sample following UV luminosity functions and measure the Lyα emission with MUSE using the latest (second) data release from the MUSE Hubble Ultra Deep Field Survey. 
  
      
   We derive the redshift evolution of X_ LAE for M_1500∈ [-21.75;-17.75] for the first time with a equivalent width range EW( Lyα)≥65 Å and find low values of X_ LAE≲30% at z≲6. The best fit linear relation is X_ LAE = 0.07^+0.06_-0.03z -0.22^+0.12_-0.24. For M_1500∈ [-20.25;-18.75] and EW( Lyα)≥25 Å, our X_ LAE values are consistent with those in the literature within 1σ at z≲5, but our median values are systematically lower than reported values over the whole redshift range. In addition, we do not find a significant dependence of X_ LAE on M_1500 for EW( Lyα)≥50 Å at z≈3–4, in contrast with previous work. The differences in X_ LAE mainly arise from selection biases for Lyman Break Galaxies (LBGs) in the literature: UV-faint LBGs are more easily selected if they have strong Lyα emission, hence X_ LAE is biased towards higher values when those samples are used.       
  
   Our results suggest either a lower increase of X_ LAE towards z≈6 than previously suggested, or even a turnover of X_ LAE at z≈5.5, which may be the signature of a late or patchy reionization process. We compared our results with predictions from a cosmological galaxy evolution model. We find that a model with a bursty star formation (SF) can reproduce our observed LAE fractions much better than models where SF is a smooth function of time. 
  
   
   The MUSE Hubble Ultra Deep Field Survey
    Haruka Kusakabe<ref>e-mail: haruka.kusakabe@unige.ch
   Jérémy Blaizot<ref>
   Thibault Garel<ref>,<ref>
   Anne Verhamme<ref>,<ref>
   Roland Bacon<ref>
   Johan Richard<ref>
   Takuya Hashimoto<ref>  
   Hanae Inami<ref>
   Simon Conseil<ref>,<ref>
   Bruno Guiderdoni<ref>
   Alyssa B. Drake<ref> 
   Edmund Christian Herenz <ref> 
   Joop Schaye <ref> 
   Pascal Oesch <ref> 
   Jorryt Matthee<ref>
   Raffaella Anna Marino<ref>
   Kasper Borello Schmidt <ref>
   Roser Pelló <ref>
   Michael Maseda <ref>
   Floriane Leclercq <ref>,<ref>
   Josephine Kerutt <ref>, <ref>   
   Guillaume Mahler <ref>
   
    Received December 18, 2019; accepted March 26, 2020
====================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================






§ INTRODUCTION


In the early Universe, the first objects formed and filled the Universe with light. They ionized the neutral gas in the intergalactic medium (IGM) via a phenomenon called “cosmic reionization”. One of the candidates for the main source of reionization is star-forming galaxies, whose ionizing radiation, called “Lyman Continuum'’ (LyC, λ<912 Å), emitted from massive stars, is expected to leak into the IGM <cit.>. Another candidate is active galactic nuclei <cit.>. However they have recently been reported to contribute less than ≈10% of the ionizing photons needed to keep the IGM ionized <cit.>. Previous studies using the Gunn-Peterson absorption trough seen in quasar spectra <cit.> and in gamma-ray burst spectra <cit.> suggest that cosmic reionization was completed by z≈6. 
The Thomson optical depth of the cosmic microwave background measured by Planck suggests that the midpoint redshift of reionization (i.e. when half the IGM had been reionized) is at z≈7.7±0.7 <cit.>. 


Lyα emission is intrinsically the strongest UV spectral feature of young star forming galaxies, and galaxies with mostly detectable Lyα emission or with Lyα equivalent widths higher than ≈25 Å  are called “Lyα emitters (LAEs)”. Lyα emission is  scattered by neutral hydrogen gas () in the IGM, and, therefore, the detectability of LAEs is affected by the  gas fraction in the IGM. The redshift evolution of Lyα luminosity functions has thus been used to investigate the history of the neutral hydrogen gas fraction of the IGM <cit.>. Lyα luminosity functions can be used to compute the evolution of the Lyα luminosity density, and its rapid decline at z≳5.7 compared with that of the cosmic star formation rate density derived from UV luminosity functions is interpreted to be caused by IGM absorption <cit.>.  



Similarly, the fraction of LAEs among UV selected galaxies, X_ LAE, can also be used to probe the evolution of the  gas fraction of the IGM  <cit.>. X_ LAE has been reported to increase from z≈3 to 6 and then to drop at z>6. This has again been interpreted as a signature of the IGM becoming more neutral at z>6 <cit.>.
The LAE fraction is complementary to the test of Lyα luminosity functions (LFs) and has some advantages:
efficient spectroscopic observations as a follow-up of continuum-selected galaxies, which is insensitive to the declining number density of star forming galaxies, and rich information obtained from the spectra such as spectroscopic redshifts and kinematics of the interstellar medium <cit.>. It also enables us to solve the degeneracy between the Lyα escape fraction among star forming galaxies with different UV magnitudes and the comparison between luminosity densities of Lyα emission and UV continuum, which are obtained from the integration of UV and Lyα LFs. In addition, recently, <cit.> suggested that the UV magnitude-dependent evolution of the LAE fraction combined with the Lyα luminosity function can be used to constrain the ionization topology of the IGM and the history of reionization. 

Using X_ LAE to set quantitative constraints on the evolution of the neutral content of the IGM remains challenging. In particular, we need to understand whether observed variations of X_ LAE are exclusively due to variations in the IGM properties, or whether they can be attributed to galaxy evolution. Following the Lyα spectroscopic observations of Lyman break galaxies (LBGs) at z≈3 by <cit.> and <cit.>, <cit.> have found that X_ LAE among LBGs evolves with redshift and depends on the rest-frame Lyα equivalent width (EW( Lyα)) cut. They also show that X_ LAE depends on the absolute rest-frame UV magnitude (M_1500), so that UV-faint galaxies are more likely to show Lyα than UV-bright galaxies <cit.>.  One conclusion from these studies is that the evolution of X_ LAE with redshift is more prominent for UV-faint galaxies and low EW( Lyα) cuts.

However, several recent studies show lower values of X_ LAE for UV-faint galaxies (-20.25< M_1500<-18.75 mag) than those in the pioneering work of <cit.>. At z≈4 and z≈5, <cit.> show more than 1σ lower X_ LAE for the faint M_1500 and low EW cut (25 Å), though their result at z≈6 is consistent with that in <cit.>. <cit.> also investigate X_ LAE for UV-faint galaxies with a low EW cut, at z≈6. They obtain a low median value of X_ LAE, which is even slightly lower than the value previously found at z≈5, though their X_ LAE is consistent within 1σ. They conclude that the drop at z>6 is less dramatic than previously found <cit.>. <cit.> and <cit.> also suggest the possibility that the effect of an increase of the  gas fraction in the IGM is observed between 5 < z < 6. This would be consistent with a later and more inhomogeneous reionization process than previously thought, as has also been recently suggested by observations and simulations of fluctuations in Lyα forest <cit.>.
The parent LBG sample in <cit.> is selected with an additional UV magnitude cut on a normal LBG selection, while the parent sample in <cit.> is mostly based on photometric redshift (photo-z) even though it is regarded as an LBG sample in their paper. Therefore, the results of X_ LAE for the faint M_1500 are not yet conclusive, possibly due to different parent sample selections. Moreover, for the UV-bright galaxies, the redshift evolution of X_ LAE for a 25 Å EW cut has not been confirmed yet  <cit.>. <cit.> and <cit.> suggest that some previous results are affected by an LBG selection bias. As strong Lyα emission affects the red band, strong LAEs can be selected more easily compared to galaxies without Lyα emission at faint UV magnitudes. It results in a high LAE fraction of LBGs <cit.>. <cit.> assess UV completeness of their parent sample using UV luminosity functions and find that their 90% completeness magnitude is ≈-20 and -21 mag at z≈4 and z≈5, respectively. 



To summarize, it is important to obtain a firm conclusion about the evolution of X_ LAE in the post-reionization epoch in order to quantify the drop of X_ LAE at z > 6 and to assess the reliability of using X_ LAE as a good probe of reionization. However, although there are a number of observational studies of X_ LAE, uncertainties in the measurement and interpretation of X_ LAE are still a matter of debate <cit.>. One of the biggest problems is the LBG selection bias due to the different depths of selected bands in previous studies. It is worth pointing out that none of the previous studies were based on complete parent samples of UV faint galaxies (-20.25< M_1500<-18.75 mag). Completeness in terms of UV magnitudes, as well as homogeneously selected samples over a wide redshift range are essential for the determination of X_ LAE. In addition, we also need deep and homogeneous spectroscopic observations of Lyα emission over a wide redshift range. 

In this study, we use a combination of deep Hubble Space Telescope (HST) and Very Large Telescope (VLT)/Multi-Unit Spectroscopic Explorer (MUSE) data, to overcome these limitations and improve our knowledge of the evolution with redshift of X_ LAE among a homogeneous parent sample of UV faint galaxies. We use HST bands that are not contaminated by Lyα emission to measure UV magnitudes to avoid a selection bias. Deep and homogeneous spectroscopic Lyα observations at a wide redshift range have been achieved in the Hubble Ultra Deep field (HUDF) by VLT/MUSE <cit.> in the guaranteed time observations (GTO), MUSE-HUDF survey <cit.>. The LAE fraction has already been investigated with MUSE and HST data, using the MUSE-Wide GTO survey in the Cosmic Assembly Near-infrared Deep Extragalactic Legacy Survey (CANDELS) Deep region in <cit.>. Their sample are constructed with an apparent magnitude cut of F775W ≤26.5 mag for an HST catalog, which is roughly converted to M_ 1500≈-19–-20 mag at z≈3–6. However, the MUSE HUDF data enable us to measure faint Lyα emission even for faint UV sources (-17.75 mag) in existing HUDF photometric catalogs.

The paper is organized as follows. In Sect. 2, we describe the data, methods, and samples: our UV-selected samples, the MUSE data, the detection and measurement of Lyα emission, the calculation of the LAE fraction, and its uncertainties. Sect. 3 presents the LAE fraction as a function of redshift and UV magnitude. In Sect. 4, we discuss our results: the differences in LAE fraction from previous results, a comparison with predictions from a model of galaxy formation, and implications for reionization. Finally, the summary and conclusions are given in Sect. 5. 

Throughout this paper, we assume a flat cosmological model with a matter density of Ω_ m = 0.3, a cosmological constant of Ω_Λ = 0.7, and a Hubble constant of H_0 = 70 km s^-1 Mpc^-1 (h_100 = 0.7). Magnitudes are given in the AB system <cit.>. 





§ DATA, METHODS, AND SAMPLES

In Sect. <ref>, we first discuss the construction of a volume-limited UV-selected sample of galaxies from the HUDF catalog of <cit.>[<https://archive.stsci.edu/hlsp/uvudf>]. In Sect. <ref> we explain how we use MUSE data to detect and measure Lyα emission from galaxies of our UV sample. In Sect. <ref> we lay out our calculations of X_ LAE and discuss the error budget. In Sect. <ref> we present our measurement of slopes of X_ LAE as a function of z and M_1500. We discuss briefly in Sect. <ref> the effects of extended Lyα emission. 



 §.§ UV-selected samples


We build our parent sample of high-redshift UV-selected galaxies using the latest HUDF catalog from R15. Sources in this catalog are detected in the average-stacked image of eight HST bands: four optical bands from ACS/WFC (F435W, F606W, F775W, and F850LP), and four near infrared (NIR) bands from WFC3/IR (F105W, F125W, F140W, and F160W). In total, out of the 9969 sources in R15's catalog,1095 and 7904 objects are within the footprints of the  and  regions of the MUSE HUDF Survey[The survey consists of two layers of different depths: a shallower area with 9 MUSE pointings () and a deeper area with 1 MUSE pointing () within the . More details are given in Sect. <ref>.] <cit.>, respectively (the duplicated region in  is removed). We note that the F140W image only covers 6.8 arcmin^2 of the 11.4 arcmin^2 footprint of the R15's catalog. The F140W photometry is only used when it is available in R15. As discussed in footnote 2 of <cit.>, the lack of F140W may affect the detection of sources in R15's catalog. Moreover, the footprint of F140W is also covered by deeper NIR images (F105W, F125W, and F160W). Indeed, the fraction of sources identified by R15 at z≳6 within the footprint of the F140W image is higher than where there is no F140W data.

In order to avoid contamination from neighboring objects, we follow <cit.> and discard all HST sources which have at least one neighbour within 06. Such associations cannot be resolved in our MUSE observations where the full width at half maximum (FWHM) of the average seeing of DR1 data is ≈06 at 7750 Å. This procedure excludes ≈ 20% of the sources. We assume that this does not result in a significant bias as it is effectively only decreasing the survey area. This assumption is true if interacting systems are not more often LAEs than isolated systems.  

R15 provide photometric redshifts and associated errors for all objects. These are obtained via spectral energy distribution (SED) fitting of photometric data in 11 HST bands using either the Bayesian Photometric Redshift (BPZ) algorithm <cit.> or the EAZY software <cit.>. In the present paper, we choose to use the results from BPZ because they are found to be more accurate <cit.>. 
 Note that R15 do not include Spitzer/IRAC data in their SED fitting. We state in Appendix <ref> that this addition would not improve the photometric redshifts of the faint galaxies studied here. Below, z_ p denotes the photometric redshift given in R15[The photometric redshift value is the one which maximizes the likelihood estimate in BPZ.], and we use it to define redshift selections of all our UV samples (see Table <ref>). In R15's catalog, the 95% lower and upper limits of z_ p are provided as uncertainties on z_ p. We use these to construct an inclusive parent sample, that we will use for Lyα searches. This sample includes all sources with photometric redshift estimates (95% confidence interval) overlapping with the redshift range (2.91– 6.12) where Lyα can be observed with MUSE. Note that we remove sources at z_ > 6.12 from our sample because the parent photometric-redshift sample may be affected by selection bias (see Sect. <ref>), and because Lyα detectability in MUSE spectra is strongly reduced by sky lines <cit.>. In the end, this inclusive parent sample consists of 3233 and 402 sources in the  and , respectively (without duplication). 

 We derive the absolute rest-frame UV magnitude using two or three HST photometric points to fit a power-law to the UV continuum. The power-law describes the spectral flux density f_ν as  f_ν= f_0 (λ/λ_0)^β+2, where λ_0 = 1500 Å, f_0 is the spectral flux density at 1500 Å (in erg s^-1 cm^-2 Hz^-1), and β is the continuum slope. We then simply have M_ 1500=-2.5log(f_0)-48.6 -5log(d_ L/10) + 2.5log(1 + z_ p), where d_ L (pc) is the luminosity distance. We choose the HST filters following <cit.> so that Lyα emission and IGM absorption are not included in the photometry: we use HST/F775W, F850LP, and F105W for objects at 2.91≤z_p≤4.44; F105W, F125W, and F140W for objects at 4.44< z_p≤5.58; and F125W, F140W, and F160W for objects at 5.58< z_p≤ 6.12. While <cit.> use the MUSE spectroscopic redshifts, we use  the values of the photometric redshifts from <cit.>. Our derived M_ 1500 values are consistent with those of <cit.> for the sources which we have in common (LAEs). The standard deviation of the relative difference in M_ 1500 for the sources included in both studies is ≈3% without a systematic offset.
 
 
 In Figure <ref>, we show the distribution of M_ 1500 as a function of z_ p for sources in our parent sample that have 2.91<z_ p≤6.12. In order to construct a complete parent sample in terms of M_ 1500, we define a limiting magnitude M_1500^ lim so that objects brighter than M_1500^ lim are detected with a signal-to-noise ratio (S/N) larger than two in at least two HST bands among the rest-frame UV HST bands. To compute M_1500^ lim, we again use a power-law continuum model, this time with a fixed UV slope β=-2, which is commonly used as a fixed value <cit.>. At each redshift, we derive the normalisation M_ 1500^ lim such that the flux can be detected in at least two HST UV bands with S/N>2. The resulting M_ 1500^ lim is shown with the thick black curve over-plotted to the data in magenta in the upper panel of Figure  <ref>. From this panel, we see that we can
build a complete UV-selected sample at redshifts z≈3–4 down to M_ 1500^ lim≈-16.3 mag, and even at z≈6 we can achieve completeness down to M_ 1500^ lim≈-17.7 mag. 

In the upper panel of Figure <ref>, we highlight two regions of parameter space that we select to do two complementary studies: X_ LAE vs. M_ 1500 with z in the range [2.91;4.44] (the polygon marked with a solid black line) and X_ LAE vs. z with z in the range [2.91;6.12] (the dashed-line rectangle). To define the criteria for X_ LAE vs. M_ 1500 plots, completeness simulations for Lyα emission (see Sect. <ref>) and z_ p bins for the M_1500^ lim calculation are also taken into account. In the lower panel of Figure <ref>, we show for comparison the locus of previous studies in the (M_ 1500-z) plane. The faint galaxies in <cit.> are shown by the light-grey shaded area. <cit.> use a UV magnitude cut of F160W ≤27.5 mag at z≈6 in their sample, which is shown with a dark-grey arrow. The recent sample of <cit.> is shown with dark-grey crosses which indicate the UV magnitudes at which they reach ≈90% completeness. The LAE fraction with the MUSE-Wide GTO survey data <cit.> adopt an apparent magnitude cut of F775W ≤26.5 mag for an HST catalog <cit.>, which we roughly convert to M_ 1500 for illustrative purposes, and show with the solid grey line in the lower panel of Figure <ref>.  Our MUSE-Deep data combined with the HST catalog from R15 allow us to probe deeper than all previous work, and to extend our study to UV-faint galaxies (i.e., M_ 1500≥-18.75 mag).

In Figure <ref>, we demonstrate the completeness of our UV-selected sample at different redshifts by comparing our UV number counts to what we would expect from the UV luminosity function (UVLF) of <cit.>. We find that the distribution of M_ 1500 for our samples (magenta) follow well those expected from the UVLFs for the same survey area at similar redshifts (black solid lines) within 1σ error bars. For comparison, we also show in the bottom panels of Figure <ref> the distribution of magnitudes of the sample of <cit.>. Clearly, their samples are incomplete even at relatively bright magnitudes (≈ -20.25 mag), most likely because of the shallow depth of their data and the LBG selection. Note that the LBG samples in <cit.> consist of two LBG samples. One is from their own sample <cit.>, and the other is a sample from the literature, which is biased towards bright objects including few B and V dropouts with magnitudes fainter than 26 mag in F850LP (z) according to <cit.>. Generally, the M_ 1500 ranges used in LAE fraction studies in the literature are close to those in <cit.>. Recently, <cit.> tested the UVLF of their LBG sample (mainly constructed from photo-z samples) used in their LAE fraction study. At z≈4 and z≈5, their LBG samples follow the UVLF in <cit.> at M_ 1500≲-20.5 mag (dark-grey dashed lines in Figure <ref>). However, it is still not complete in the UV for the faint sample (-20.25≤ M_ 1500≤-18.75 mag) just as <cit.>.These comparisons illustrate the methodological improvement of our study in terms of the M_ 1500 completeness of the sample of galaxies for which we estimate the LAE fraction.

For future use, N_ 1500(z_ p,  M_1500) denotes the number of galaxies with a photometric redshift z_ p and absolute magnitude M_1500 within the given ranges (see Table <ref> which summarises different samples). As discussed above, our UV-selected samples are volume-limited and N_ 1500 is directly measured from the catalog with no need for incompleteness corrections. 

 


 






 §.§ Counting LAEs within the UV-selected sample
 


  §.§.§ MUSE data


The data of the MUSE HUDF Survey were obtained as part of the MUSE GTO program (PI: R. Bacon). The MUSE HUDF Survey design is presented in <cit.>. It consists of two layers of different depths: the  is composed of 9 MUSE pointings that cover a 3'× 3' area (9.92 arcmin^2) with an integration time of ≈10 hours; the  is a deeper integration at a 1'× 1' sub field within the , with an integration time of ≈31 hours. MUSE covers a wide optical wavelength range, from 4750 Å to 9300 Å, which allows the observation of the  line from z≈2.9 to z≈6.6. The typical spectral resolving power is R=3000, with a spectral sampling of 1.25 Å. The spatial resolution (pixel size) is 02 × 02 per pixel.

In the present paper, we use the latest data release (the second data release, hereafter DR2) from the MUSE HUDF (Bacon et al.  in prep.). The improved data reduction process results in data cubes with fewer systematics and a better sky subtraction. The FWHM of the Moffat point spread function (PSF) is 065 at 7000 Å in the MUSE HUDF. The estimated 1 σ surface brightness limits are 2.8×10^-20 erg s^-1 cm^-2 Å^-1 arcsec^-2 and 5.5×10^-20 erg s^-1 cm^-2 Å^-1 arcsec^-2 in  and , respectively, in the wavelength range of 7000–8500 Å (excluding regions of OH sky emission, see , Bacon et al. in prep. for more details).  For instance, the estimated 3 σ flux limits are 1.5×10^-19 erg s^-1 cm^-2 and 3.1×10^-19 erg s^-1 cm^-2 in  and , respectively, for a point-like source extracted over three spectral channels (i.e., 3.75Å) around 7000 Å <cit.>. The PSF and noise characteristics are similar to the DR1 data, except in the reddest part of the wavelength range.

In order to measure the fraction of galaxies which have a strong  line, we first extract a 1D spectrum from the MUSE cube for each HST source in our parent sample. We proceed as follows. First, we convolve the HST segmentation map of the R15 catalog with the MUSE PSF, which is normalized to 1. To obtain a spatial mask applicable to MUSE observations for each object, we apply a threshold value of 0.2 to the normalized convolved segmentation map. The median value of the radius of the normalized mask is ≈10 to 08 arcsecond at z≈3 to 6, which is not affected by Lyα halo flux <cit.>. Second, we integrate the cube spatially over the extent of the mask. We note that PSF weighted or white-light weighted integrations are used to extract spectra in the DR1 catalog of <cit.>. These provide a higher S/N in the extracted spectrum. However, in the present paper, we do not use a spatial weighting. This results in slightly lower S/N values, but more accurate estimates of the fluxes (i.e., conserved flux), which are needed to assess the completeness of our  detections. Third, we subtract local residual background emission from the extracted spectrum for the 1D spectra as in <cit.>. The local background is defined in 5 × 5 subcubes avoiding the masks of any source.

 


  §.§.§ Search for Lyα emission



In order to detect Lyα emission lines in the 1D spectra extracted above, we use a customized version of the  software[The original  in <cit.> is based on a cross-correlation algorithm <cit.> and is publicly available at: <https://github.com/Samreay/Marz.>] <cit.> described in <cit.>.  compares 1D spectra to a list of templates and returns the best-fitting spectroscopic redshift, the best-fitting 1D template, and a confidence level for the result (called the quality operator, QOP). In our customized  version, the list of templates consists of templates made using MUSE data, and the interface is improved <cit.>. We use our version of  in a similar manner to <cit.>, except for the two following changes. First, we do not activate cosmic ray replacement in  because (1) it affects the detectability of bright and spectrally peaky Lyα emission, and (2) cosmic rays are efficiently removed in the data reduction. Second, we only use template spectra with Lyα emission: those of IDs=10, 18, and 19, which are used in <cit.>, and those of IDs=25, 26, 27 and 30, which are newly built from MUSE data in Bacon et al. in prep. and show single peaked Lyα (see Appendix <ref> for the template spectra). 
As in <cit.>, we use the 1D spectra and source files including the subcubes and cutouts of HST UV to NIR images for the parent sources as input for . 


To select robust Lyα detections, we keep only galaxies which  identifies as LAEs with a high confidence level <cit.>[The confidence level is given as QOP, which is calculated from the peak values of the cross-correlation function <cit.>. In our version of  <cit.>, QOP =3 (QOP =2) is regarded as “Great” (“Good”) and  corresponds to 99.55% (95%) confidence in the original  <cit.>. However, the original relation between QOP and confidence percentage is calibrated with SED templates different from ours, and the confidence percentage may not be directly applicable to our data. Note that the FOM criterion for QOP =3 (QOP =2) in our version of  is the same as that for QOP ≥4 (QOP =3) in the original  (see Figure1 in  and Figure 12 in ). ]. Sources with lower confidence levels are not regarded as detected LAEs. According to this selection, among the 3233 (402) sources in   (), 374 (70) are LAEs. However, some of these LAE candidates are in fact  emitters or non-LAEs polluted by extended Lyα emission from LAE neighbors. We visually inspect all the LAE candidates as in <cit.>: we check the entire MUSE spectra, Lyα line profiles, MUSE white-light images, MUSE narrow-band images of Lyα emission, all the existing HST UV to NIR images, and HST colors by eye using the customized . The MUSE white-light image is created by collapsing the 5 × 5 MUSE subcubes in wavelength direction <cit.>, while the narrow-band image for the Lyα emission is extracted from the wavelength range around the Lyα emission in the subcubes for the sources <cit.>. In contrast to <cit.>, we also use a consistent photometric redshift (95% uncertainty range) as an evidence of Lyα emission. As a result, we have 276 (58) LAEs at z≈2.9–6.1 among the parent sample in the  () field. Most of the removed sources (≈80%) have a 1D spectrum contaminated by (extended) Lyα emission from neighboring objects, which can be distinguished using the Lyα narrow-band images, MUSE white-light images, and HST images. We show an example of Lyα contamination in Appendix <ref>. 





  §.§.§ Measurement of Lyα fluxes
 

For our LAEs, we measure the Lyα fluxes from the 1D spectra used for Lyα detection and described in Sect. <ref>. The aperture size is defined by the R15's segmentation map for each source convolved with the MUSE PSF (see Sect. <ref> and Sect. <ref> for more details). 
It has been reported that Lyα emission is often spatially offset from the stellar UV continuum <cit.>. The typical offset values for our LAEs are however measured to be less than 01 <cit.>, significantly less than the PSF scale of our observations. We therefore assume that all the flux of the central Lyα component can be captured by our apertures, centered on the continuum emission peak, and with median radius ranging from ≈10 to 08 at z≈3 to 6. We fit the Lyα emission either with an asymmetric Gaussian or a double Gaussian profile. The choice between these two solutions is made after visual inspection. In practice, we use the  and  methods of the publicly available software  <cit.>[MPDAF, the MUSE Python Data Analysis Framework, is publicly available from the following link: <https://mpdaf.readthedocs.io/en/latest/>]. We use the spectroscopic redshift from  as information for the center wavelength of the fit, with a fitting range of rest-frame 1190 Å to 1270 Å. We inspect all the Lyα spectral profile fits to confirm their validity.

Once we have measured the Lyα fluxes, we compute the rest-frame Lyα equivalent width using M_ 1500 and β defined in Sect. <ref> to estimate the continuum. The distribution of EW( Lyα) as a function of redshift is shown in Figure <ref>.


 





  §.§.§ Completeness estimate and correction
 

 
 
In order to estimate the detection completeness of Lyα emission for the MUSE HUDF data with , we insert fake Lyα emission lines into 1D spectra and try to detect them as explained in Sect. <ref>. 

We take realistic noise into account by creating 1D sky-background spectra from MUSE sub-cubes extracted for different continuum-selected sources as detailed in Sect. <ref>. We choose a sample of spectra which show a clear Lyα emission line <cit.>, and which do not have continuum or other spectroscopic features. We mask the spectral pixels covered by the Lyα emission in these 1D spectra (including ± 20 pixels ≈± 25 Å  around the line center). Note that we do not insert fake Lyα lines in the masked regions. With this procedure, we obtain 131 (35) 1D sky-background spectra in the  () field.

We add fake emission lines with fluxes taking 18 values regularly spaced in log between 6 × 10^-19 erg s^-1 cm^-2 and 50 × 10^-18 erg s^-1 cm^-2, and 3102 redshift values regularly distributed between z=2.93 and z=6.12. For each flux-redshift pair, we draw 4 lines (yielding a total of 223,344 lines) and add each of them to one of our 131 (35) 1D spectra chosen randomly. Each fake Lyα line, has line-shape parameter values (total FWHM and FWHM ratio of red wing to blue wing) randomly drawn from the measured distribution of Lyα emission line shapes of LAEs used in <cit.> and <cit.>. We use the  method of  to generate the fake lines and add them to our test spectra. 

We then repeat the detection procedure of Sect. <ref>, applying the same cut at a confidence level of “Good”. In each field ( or ), we compute the completeness of Lyα detection as a function of the Lyα flux f_ Lyα for five redshift bins: 2.91≤ z <3.68 (z≈3.3), 3.68≤ z <4.44 (z≈4.1), 4.44≤ z <5.01 (z≈4.7), 5.01≤ z <5.58 (z≈5.3), and 5.58≤ z ≤6.12 (z≈5.9), which are defined from the redshift bins used to derive M_1500. We fit each simulated completeness curve with a formula based on the error function, <cit.>:

    C(f_ Lyα)=1/2[1- erf( -2.5log_10f_ Lyα -a/√(2b))],

where a and b are two free parameters for fitting (see Figure <ref>). We use the function  from  to perform the fit. The best fit parameters for the completeness curve in the  and   are summarized in Table <ref>. 

At completeness above ≈0.8, our best fit relations slightly overestimate the measured completeness. The analytic fit is however at most ≈5% above the 1σ upper errors, and this does not have a noticeable effect on the calculation of X_ LAE. We nevertheless take this into account in the error propagation of X_ LAE in Sect. <ref>.  

Theoretically, completeness functions should just scale with S/N and thus be applicable throughout the wavelength (or redshift) range of the instrument. For the MUSE-WIDE survey, <cit.> indeed find that the shape of their completeness function is  independent of redshift. As expected, we also find very similar behavior of completeness at z≈3.3 and 4.1, in both fields. At these redshifts, the noise is well behaved and there are only few accidents due to sky-line removal in the spectra. At z≈ 4.7, the shape of the completeness curve is still well described by Eq. <ref>, but the curve is shifted to fainter flux with a shallower slope. At z≈5.3 and 5.9, the shapes of the best-fitting completeness are different from those at lower redshifts: they have a shallower slope, the data points with completeness above 0.8 are not well fit, and the completeness at a given flux is much lower than that at lower redshifts. The lower normalisation and distorted shape of the completeness may be caused by the many sky emission lines at high redshifts <cit.>. Because  is not a local line detector, as opposed, e.g. to the matched-filtering approach implemented in the tool  utilized in <cit.> where the filter has a compact support in spectral space,  is affected by relatively long-range noise or distant spectroscopic features in the spectra. Thus,  often does not return a high confidence level (“Great” and “Good”) for LAEs at z≳5, and our completeness goes down at z≳5 even for relatively bright Lyα fluxes. We note that the LAE template spectra that we use all have a single-peaked Lyα profile (see Figure <ref>). However, the exact shape of the line profile is shown to have little impact on the detectability with cross-correlation function in general <cit.>. The shape of the Lyα line of 's template should not affect significantly the detection rate of LAEs (see also Appendix <ref>). Indeed, for example, our sample contains LAEs with double-peaked lines even though none of our templates have such features. We experimented using all the templates in , including templates of different galaxy populations such as  and  emitters, and we found only a small impact on our LAE sample, which is well within the error bars. Finally, we checked the dependence of completeness on the FWHM of fake Lyα emission lines at fixed flux, and again found no significant trend. 







 
 
 


In the following, N_ LAE^ det(z_ s, M_1500,  EW) denotes the number of galaxies with detected Lyα emission and with spectroscopic redshift z_ s, absolute magnitude M_1500, and rest-frame equivalent width EW in given ranges. We estimate the true number of LAEs with a corrected value N_ LAE^ corr defined as follows. For a given field and redshift bin, we use the fits to the completeness function above to define four Lyα flux bins which correspond to regularly spaced bins in the logarithm of the completeness (C), ranging from C=0.1 to C=0.9. We then count the number of detected LAEs (within a given z_ s, M_1500, and EW bin) in each flux bin and divide it by the mean completeness (in log) in each of the flux bins. We compute N_ LAE^ corr as the sum of these over the four flux bins. When the uncertainties of the LAE fraction is calculated, the completeness correction value in each flux bin is propagated, and the uncertainty of completeness correction itself is taken into account as described in the next section. The number of flux bins is defined through a test described in Appendix <ref>. Four to six bins are a sweet spot where the error bars are small and they appear converged. For a larger number of bins, we often get flux-bins with no object at all, and these bins contribute to a large error bar. For a smaller number of bins, we introduce a larger error on the completeness correction (averaged over the bin). We adopt four bins. 




 §.§ X_ LAE and its error budget
 

Knowing the number of UV-selected galaxies in a volume-limited sample (N_ 1500(z_ p,  M_1500), Sect.<ref>), and the number of LAE among the inclusive parent sample (N_ LAE^ corr(z_ s, M_1500,  EW), Sect. <ref>), the fraction of UV-selected galaxies with a Lyα line is simply given as:

    X_ LAE = N_ LAE^ corr(z_ s, M_1500,  EW)/N_ 1500(z_ p,  M_1500).


The uncertainties on X_ LAE arise from four components: (1) the uncertainty due to contaminants (type II error) and missed objects (type I error) in N_ 1500(z_ p,  M_1500),  (2) the uncertainty due to the completeness correction of N_ LAE^ corr(z_ s, M_1500,  EW), (3) the uncertainty for Bernoulli trials (i.e., fraction of N_ LAE^ corr(z_ s, M_1500,  EW) over N_ 1500(z_ p,  M_1500)) measured by a binomial proportion confidence interval and (4) the uncertainty due to cosmic variance. We note that there is no obvious sample selection bias in our UV galaxies as shown in Figure <ref>, and our Lyα measurements are  homogeneous over the whole redshift range (see discussion in Sect. <ref>). We estimate the relative uncertainty of X_ LAE from error propagation and discuss each of these contributions below.  




  (1) Contaminants and missed objects in N_ 1500(z_ p,  M_1500)
 
Sources with z_ p in a given z range that are truly located outside of the z range are contaminants in the parent continuum-selected sample, while sources with z_ p outside of the z that are truly located in the z range are missed sources. These mismatches of z_ p can happen because of confusion between Lyman and 4000 Å breaks in the SED fitting. In addition, IGM absorption modeling has been suggested to affect z_ p estimation <cit.>. 
As discussed in <cit.> and <cit.>, the fraction of contaminants is very low for high confidence level objects <cit.>. The fraction of missed objects with a relative redshift difference of more than 15%, |z_ s-z_ p|/(1+z_ s)>0.15, is suggested to be ≈10% <cit.>. Since the missed objects whose 95% uncertainty range for z_ p is outside the z range for MUSE-LAEs are not included in our parent continuum-selected sample, they are also not included in our LAE sample. With an assumption that the fractions of missed objects are the same for the parent and LAE samples, the uncertainties due to missed objects can be neglected as well as those due to contaminants. Note that we do not find any significant relation between the z_ p-z_ s difference and Lyα EW as well as Lyα flux. 


  (2) Completeness correction of N_ LAE^ corr(z_ s, M_1500,  EW)
 
Simulated data of completeness are fitted well with Equation (<ref>) at z≈3.3, 4.1, and 4.7. Even at z≈5.3 and 5.9, the differences in completeness between the simulated data and the best-fitting functions are at most ≈5%, which is much smaller than the uncertainty due to the flux binning. The flux binning for the completeness correction described in Sect. <ref> causes an uncertainty of at most ≈±32%. The completeness bins, spaced regularly in log from 0.1 to 0.9, correspond to steps of a factor 1.73, which has a square-root of ≈ 1.32. 
We use this very conservative estimate of 32% for the completeness correction error in our error budget. The completeness correction error is smaller than the error component (3) as described later and does not change the total uncertainty of X_ LAE. We note that the completeness correction value is also taken into account according to error propagation, when the error component (3) is calculated.



   (3) Uncertainties for Bernoulli trials
 
Measuring the fraction of a sub-sample among a parent sample is a kind of experiment of Bernoulli trials. An uncertainty for a Bernoulli trial is given by a binomial proportion confidence interval (hereafter, BPCI). We use the python module  from  and provide an approximate uncertainty for a given confidence interval (=68%, 1σ, in this work), the number of trials, and the number of successes. Here, the number of trials and the number of successful experiments are N_ 1500(z_ p,  M_1500) and N_ LAE^ corr(z_ s, M_1500,  EW), respectively. However, we cannot obtain N_ LAE^ corr(z_ s, M_1500,  EW) directly from the observations. To include the effect of the completeness correction in each Lyα flux bin (described in Sect. <ref>) in the error propagation, we calculate the uncertainty of the LAE fraction in each Lyα flux bin without applying a completeness correction from  and multiply the uncertainties by the correction value in the flux bin. We choose as an approximation formula the Wilson score interval <cit.>, which is known to return an appropriate output even for a small number of trials and/or success experiments. Our method is confirmed to be accurate by numerical tests described in Appendix <ref>. For the flux-bins with no LAEs, the average completeness value among the bins and the number of LAEs (=0) are used to derive the uncertainties conservatively. When we sum over all the uncertainties in a flux bin to derive the total uncertainties for the component (3), a python module,  developed in <cit.>, is used to treat asymmetric errors by BPCI. We note that the Poisson errors of N_ LAE^ corr(z_ s, M_1500,  EW) and N_ 1500(z_ p,  M_1500) are commonly used in the literature. However, the error for the LAE fraction should be derived by BPCI like in other fraction studies <cit.> to obtain statistically correct errors. The BPCI method is reviewed for astronomical uses in <cit.>.



  (4) Cosmic variance
 
The survey volume in each redshift range is limited to ≈ 1.5-2.5 × 10^4 cMpc^3. However, we find that the uncertainty due to cosmic variance is less than the BPCI error and thus not affecting our X_ LAE significantly (See Sect. <ref> for more details). Since the uncertainty due to cosmic variance  cannot be included in our MUSE measurements, we neglect this error component (4). 


In addition to (1) - (4), uncertainties in photo-z estimations and flux measurements are also potential error components. In Appendix <ref>, we discuss the impact of uncertainties on z_ p and conclude that it does not affect the error bar of X_ LAE significantly. Uncertainties of M_1500, β, and the Lyα flux, combine into uncertainties on EW( Lyα). As shown in Figures <ref> and <ref>, X_ LAE shows a slight dependence on M_1500 and EW( Lyα). We thus expect that a small error on these quantities will translate into an even smaller error on X_ LAE.  Note that although some objects have large errors in M_1500 in Figure <ref>, they are not included in our analysis because they are very faint[Among the subsamples shown in Table <ref>, fainter-M_1500 and higher-z_ p objects have greater uncertainties in M_1500. The medians (standard deviations) of the uncertainties for the subsamples with -18.75≤ M_1500≤-17.75 mag are 0.05 mag (0.02 mag) at z≈3.3 and 0.15 mag (0.12 mag) at z≈5.6. Those for the subsamples with -19.0≤ M_1500≤-18.0 mag is 0.04 mag (0.01 mag) at z≈3.3 and 0.08 mag (0.11 mag) at z≈4.1. These uncertainties are much smaller than the width of M_1500 bins. With regard to S/N cuts of Lyα fluxes corresponding to EW( Lyα) cuts, those in the  field for M_1500=-17.75 mag and EW( Lyα)=65 Å are estimated to be ≈17.4 at z≈3.3 and ≈4.4 at z≈5.6, if we assume β=-2. Similarly, the S/N cuts for M_1500=-18.0 mag and EW( Lyα)=50 Å are ≈16.9 at z≈3.3 and ≈11.1 at z≈4.1. Here, the noise is the median of those shown in Figure <ref> in each redshift bin.]. Thus, we ignore these two uncertainties, as is commonly done in the rest of the literature. 





With these considerations, the uncertainty of the LAE fraction is the quadratic sum of the uncertainty terms (2) and (3). Below, the error bars on X_ LAE represent the 68% confidence intervals around the values calculated by Equation (<ref>). Note that the dominant error for X_ LAE derived in this work is component (3), the uncertainties for a Bernoulli trial, which are, for instance, 38% and 78% of X_ LAE for -21.75≤ M_1500≤-18.75 mag and EW( Lyα)≥65 Å at z≈3 and 5.6, respectively.




 §.§ Measurement of the slope of X_ LAE as a function of z and M_1500
  
We measure linear slopes of X_ LAE as a function of z and M_1500 using a python package for orthogonal distance regression (hereafter, ODR) fitting, , to account for widths of bins in x-axis and uncertainties of X_ LAE in y-axis. The ODR fitting minimizes the sum of squared perpendicular distances from the data to the fitting model. Since uncertainties of X_ LAE are not symmetric, a Monte Carlo simulation with 10000 trials is used. We assume an asymmetric Gaussian profile as a probability distribution function for X_ LAE with each of the upper and lower uncertainties at a given bin (z or M_1500) in x-axis. We fit a linear relation (y=ax+b) in each trial drawing X_ LAE randomly with . The best fit values of a and b and their error bars are derived from the median values and the 68% confidence intervals around the median values. The results of fitting are shown in Sect. <ref>. 





 §.§ Extended Lyα emission
 
Our aim is to measure how the fraction of UV-selected galaxies showing a strong  line varies with redshift. We make the choice to discard possible significant contributions to the  luminosities by extended  haloes <cit.> in our study, though the contribution of LAHs to the total Lyα fluxes is typically more than ≈50% <cit.>. There are a number of reasons for this choice. First, it is largely motivated by our lack of understanding of the physical processes lighting up these halos. In particular, it is not clear how they relate to the UV luminosities of their associated galaxies <cit.>, to what extent they are associated to star formation <cit.>, or whether the nature of this association could vary with redshift. In order to assess the evolution of X_ LAE with redshift, it thus appears more conservative to limit our measurement of the  emission from galaxies to the part which most likely has the same origin as the continuum UV light. Any evolution is then more likely to be related to the evolution of the ionisation state of the IGM. With the above procedure, our 1D spectra include as little as possible of the extended  emission that is found around LAEs <cit.>. Second, our choice has the advantage of following a similar methodology without including halo fluxes used in other studies <cit.>, and thus allows for a fair comparison. Third, it is difficult to measure the faintest halos. If we demanded that LAHs were detected around galaxies in our sample, we would limit our sample to the brightest and/or compact halos only <cit.>. So even though an IFU enables us in principle to separate the central and halo component more clearly than slit and fiber spectrometers, the signal-to-noise ratio required to do so is still prohibitive for statistical studies such as ours. 

Note that because of our choice, the Lyα fluxes and EWs in the present paper are smaller than the total Lyα fluxes and EWs reported by e.g.  <cit.>, <cit.>, and  <cit.>.







§ RESULTS


In order to measure the variation of X_ LAE with redshift or UV absolute magnitude, we design several sub-samples shown in Table <ref>. We use EW( Lyα) cuts starting at 25 Å, which is a common limit in the literature, and then increase in steps of 20Å to 45 Å, 65 Å, and 85 Å. We also use 50 Å and 55 Å cuts, for comparison to <cit.> and <cit.>. In Sect. <ref>, we present our results for the X_ LAE-z relation, going as faint as M_1500=-17.75 mag for the first time in a homogeneous way over the redshift range z≈3 to z≈6. We discuss how these results compare to existing measurements. In Sect. <ref>, we present the first measurement of the X_ LAE-M_1500 relation for galaxies as faint as M_1500=-17.00 mag at z≈3, and compare our findings to other studies. The numerical values of X_ LAE are summarized in Tables <ref> and <ref>. The slopes of the best fit linear relations of X_ LAE as a function of z and M_1500 are shown in Figure <ref> in Appendix <ref>, and summarised in Tables <ref> and <ref>. 













 §.§ Redshift evolution of X_ LAE


We derive the redshift evolution of X_ LAE for EW( Lyα)≥65 Å from M_1500=-21.75 mag to the faint UV magnitude of -17.75 mag. We show the results in the upper panel of Figure  <ref> with the large filled hexagons. We find a weak rise of X_ LAE from z≈3 to z≈6, even though poor statistics do not allow us to set a firm constraint at z≈5.6. Breaking our sample into a bright end (with -21.75≤ M_1500≤-18.75 mag, purple circles), and a faint end (with -18.75≤ M_1500≤-17.75 mag, purple squares), we find similar trends that are consistent within the 1σ error bars, suggesting that the X_ LAE-z relation does not depend strongly on the rest-frame UV absolute magnitude. The best fit linear relations are X_ LAE = 0.07^+0.06_-0.03z -0.22^+0.12_-0.24, X_ LAE = 0.05^+0.05_-0.03z -0.13^+0.12_-0.18, and X_ LAE = 0.09^+0.09_-0.04z -0.28^+0.18_-0.36 for -21.75≤ M_1500≤-17.75 mag, -21.75≤ M_1500≤-18.75 mag, and -18.75≤ M_1500≤-17.75 mag, respectively. Note that the bright sample here is dominated by the more numerous sub-L^∗ galaxies, which are fainter than the bright samples in the literature.

If we select the brighter part of our sample (-21.75≤ M_1500≤-18.75 mag), we can estimate X_ LAE down to lower equivalent widths. In the lower panel of Figure <ref>, the orange circles show X_ LAE for galaxies with EW( Lyα)≥45 Å. This is also found to increase from z≈3 to z≈4–6, and is above the relation obtained for EW( Lyα)≥65 Å also shown in the lower panel, as expected. The best fit linear relation is X_ LAE = 0.05^+0.05_-0.03z -0.07^+0.14_-0.20. These results are qualitatively consistent with the trend of increasing X_ LAE with increasing z based on bright samples in the literature (see below). 

The fact that only ≈0–30% of galaxies within the UV magnitude range at z≤5 are observed as LAEs with EW( Lyα)≥65 Å requires some explanation of the physical mechanisms. We discuss this further in Sect. <ref>.






Next, we compare our MUSE results of X_ LAE with previous studies. For this purpose, we derive X_ LAE for M_1500 in the range [-20.25;-18.75] <cit.>, and for EW( Lyα)≥25 Å and EW( Lyα)≥55 Å. Figure  <ref> shows our results and those of other studies as a function of redshift. At z≲5, we confirm the low values from <cit.> (grey crosses at z≈4 and z≈5) for EW( Lyα)≥25 Å. Our median values of X_ LAE at z≈4.1 and z≈4.7 are somewhat smaller than those at z≈4 and z≈5 from <cit.>, although they are compatible within the large error bars. At z≈5.6, our value for EW( Lyα)≥25 Å appears to be significantly lower than those reported in the literature. We note that the discrepancy between our work and <cit.> is however only 1.14σ, and might thus be caused by the statistical error. However, our result is more than 2σ away from those of <cit.> and <cit.>, which is less likely to be statistical fluctuation. We discuss in Sect. <ref> potential biases which may explain why these two latter references find large values of X_ LAE. We discuss the effect of cosmic variance in Sect. <ref>, and find that it cannot explain such large differences. Another possibility is that the low value we find reflects a late and/or patchy reionization process, and we discuss that further in Sect. <ref>.


We also check the best fit linear relation of X_ LAE as a function of z for EW( Lyα)≥25 Å, which is X_ LAE = 0.01^+0.05_-0.05z +0.16^+0.20_-0.22. The best fit slope is lowered by the point at z≳5.6, and is shallower than 0.11±0.04 in <cit.> and 0.18^+0.06_-0.06 in <cit.>. It is consistent with the flat relation reported by <cit.> within the 1σ error bars (0.014±0.02) and that by <cit.>. Note that <cit.> discuss the slope of X_ LAE against z using a sample of the MUSE-Wide GTO survey with an apparent magnitude cut of F775<26.5 mag, which is shown by a grey solid line in Figure <ref>. They also include the contribution of extended Lyα halos to their Lyα fluxes, which enhances the values, contrary to us. 
With regard to EW( Lyα)≥55 Å  the best fit relation is X_ LAE = 0.04^+0.05_-0.03z -0.05^+0.14_-0.19, whose slope is consistent with that in <cit.>, 0.018±0.036. 









 §.§ UV magnitude dependence of X_ LAE
 

Figure  <ref> shows a diagram of X_ LAE and M_1500 at z=2.91–3.68 (≈3.3) for our MUSE sample. This is the first time that the dependence of the LAE fraction on M_1500 is studied at M_1500≥-18.5 mag. The LAE fractions for EW( Lyα)≥25 Å (45 Å, 65 Å, and 85 Å) are shown with the purple (violet, orange, and yellow) stars. The best fit linear relations are 
X_ LAE = -0.02^+0.07_-0.08M_1500 -0.30^+1.39_-1.62,
X_ LAE = 0.01^+0.02_-0.03M_1500 +0.33^+0.45_-0.54,
X_ LAE = 0.02^+0.01_-0.01M_1500 +0.41^+0.18_-0.30, and
X_ LAE = 0.01^+0.01_-0.01M_1500 +0.12^+0.13_-0.14 for EW cuts of 25 Å, 45 Å, 65 Å, and 85 Å, respectively. We find no clear dependence of X_ LAE on M_1500 in tension with the clear rise of X_ LAE to faint UV magnitude for an EW cut of 50 Åreported in <cit.>. 


Our results also show that the LAE fraction is sensitive to the equivalent width selection, as expected e.g. from <cit.>. Although this means that the LAE fraction is useful in itself to test cosmological galaxy evolution models <cit.>, it also raises concern for the usage of X_ LAE as a probe of the IGM neutral fraction at the end of reionization <cit.>, since homogeneous measurements of Lyα emission over a wide redshift range are required for a fair comparison.  








In Figure <ref>, our results for the relation between X_ LAE and M_1500 for EW( Lyα)≥50 Å at z≈3-4 (filled and open red stars) are compared with those in <cit.> (filled grey squares). The best fit linear relations for our results at z≈2.9–3.7 and at ≈3.7–4.4 are X_ LAE = 0.01^+0.02_-0.03M_1500 +0.34^+0.44_-0.55 and X_ LAE = 0.01^+0.06_-0.10M_1500 +0.33^+1.12_-1.80, respectively. We find no dependence of X_ LAE on M_1500 as opposed to the claim in <cit.>, whose best fit relation is X_ LAE = 0.13^+0.03_-0.03M_1500 +2.87^+0.74_-0.72. Our X_ LAE is lower than that in <cit.> at UV magnitude fainter than M_1500≈-19 mag, and possibly at M_1500≈-20 mag. We discuss the difference in X_ LAE between this work and <cit.> in Sect. <ref>.






§ DISCUSSION


In this section, we assess the cause of the differences between our results and previous results, compare our results with predictions from a cosmological galaxy formation model, and discuss the evolution of the LAE fraction and implication for reionization. 




 §.§ Possible causes of the differences between our MUSE results and previous results





In Figure <ref>, we find that our measurements of X_ LAE are systematically lower than those of <cit.>, although consistent within the error bars. This tension supports the results of <cit.>, who also find low median values of X_ LAE at z≈4–5. It is worth discussing the potential origins of this tension, since the median values of X_ LAE have been used to assess cosmic reionization in theoretical studies <cit.>. The difference between our study and that of <cit.> is more striking in Figure <ref> which shows X_ LAE as a function of M_1500. Here, our results are inconsistent with theirs at a faint UV magnitude, even when taking into account the large error bars. Below we discuss two possible origins of this discrepancy in the plot of X_ LAE as a function of M_1500: the LBG selection bias, and systematics due to different observing methods. 



  §.§.§ LBG selection bias

There is a possibility that the LBG sample of <cit.> is biased towards bright Lyα emission i.e., higher X_ LAE, as pointed out in previous studies <cit.>. In other words, LBG selections could be biased towards having higher X_ LAE if they preferentially miss low-EW sources. The LBG selection consists of a set of color-color criteria and signal-to-noise cuts <cit.>. <cit.> obtain a relatively low median X_ LAE at z≈6 and discuss the causes. They use common selection criteria for the i-dropout <cit.>, but add an additional criterion of the H-band (F160W, rest-frame UV at ≈6) magnitude cut. When Lyα emission is located in the red band, strong Lyα emission in a UV spectrum can significantly enhance the Lyman break. The additional criterion in <cit.> can suppress the LBG selection bias, which increases X_ LAE for faint UV sources as will be discussed below.

Here we estimate the effect on X_ LAE of two aspects of the LBG selection bias for B (F435W)-dropouts: the impact of Lyα contamination on the signal-to-noise ratio cut in the V (F606W) band, and on the color-color criteria in a diagram of F435W-F606W vs. F606W-F850LP (black solid line in Figure <ref>). We estimate F606W magnitudes assuming a power-law spectrum with a UV slope of -2 (λ≥912 Å in rest frame) and IGM transmission following <cit.>. As shown in the upper panel in Figure <ref>, strong Lyα emission can increase the flux in F606W noticeably, especially at z≳4. The magnitude shift becomes larger at higher redshifts because of the increasingly smaller rest-frame wavelength range of the UV continuum that is covered by F606W as the redshift goes up. Even at z≈4, however, a moderate Lyα emission line with EW( Lyα)=50 Å can cause a ≈0.16 mag difference in F606W. This affects both the signal-to-noise ratio and the colors. 

We first illustrate the effect on the signal-to-noise ratio cut of F606W by considering an object with M_1500=-19 mag. At z=4.5 (z=4.0), a source without Lyα emission has F606W=28.3 mag (F606W=27.5 mag), while a source with EW( Lyα)=50 Å has a 0.27 (0.16) brighter magnitude. Since these F606W magnitudes are close to the 5σ limiting magnitude of 28.0 mag in <cit.>, which corresponds to a completeness of 50% in the case of a S/N≥5 cut, the completeness for their B-dropout changes drastically around 28.0 mag.

Second, the effect on the two colors for the B-dropout selection is shown in the lower panel of Figure <ref>. The green, red and violet arrows show color-color shifts in the diagram of F435W-F606W and F606W-F850LP due to the magnitude shift of F606W in the case of EW( Lyα)=25 Å, 50 Å, and 100 Å, respectively. The shift of F606W enhances the possibility of meeting the dropout criteria. Indeed, at z=3.5–4.5, all of our LAEs with an absolute magnitude of F775W=-20.25–-18.75 mag are located in the dropout selection region (upper left region of bottom panel in Figure  <ref>). However, ≈10% of continuum-selected galaxies do not meet the dropout criteria. Therefore, strong Lyα emission can enhance the probability to meet LBG-selection criteria both in terms of the signal-to-noise cut and of the color-color criteria. 

We estimate the completeness of the B-dropout galaxies in <cit.> using a plot of surface number density as a function of UV magnitude for B-dropouts in <cit.>.
<cit.> use the same color-color criteria as <cit.> but with ≈0.6 mag shallower data sets than those in <cit.>.  Figure 1 in <cit.> shows the surface number density of the B-dropouts as a function of apparent F775W magnitude (i.e., apparent rest-frame UV magnitude). At absolute UV magnitudes of ≈-18.8 and -18.3 at z≈4, the completeness values are ≈25% and ≈5%, respectively. The completeness values for ≈0.6 mag shallower data in <cit.> are estimated to be ≈25% and ≈5% at M_1500≈-19.4 mag and M_1500≈-18.9 mag, respectively, if the behavior of completeness as a function of S/N is similar. As shown in the lower panels of Figure <ref>, the B-dropout galaxies in <cit.> are not complete at M_1500≈-19.0 mag.  Moreover, <cit.> adopt stricter signal-to-noise ratio cuts than <cit.>'s criteria for B-dropouts, and the completeness values in <cit.> may be lower than estimated here,  especially for faint sources. 

Following the discussions above, at z≈3.5–4.5, the observed number (the denominator of X_ LAE) for their B-dropouts with M_1500≈-19 mag is estimated to be less than 25% of the true value, while the observed number (the numerator of X_ LAE) for LAEs with EW( Lyα)=50 Å is estimated to be larger than 50% of the true value under the assumption that all the LAEs meet the color-color criteria. This means that X_ LAE for their B-dropout sample may be more than ≈1.5 times larger than the X_ LAE for a complete sample, 0.3^+0.08_-0.09. If the overestimate would be corrected for X_ LAE at M_1500≈-19 mag, their X_ LAE would be consistent with ours within the 1σ error bars as shown by the open grey square in Figure <ref>.

Therefore, the B-dropout selection bias may be the dominant cause of the difference in X_ LAE between LBGs and photo-z selected galaxies at z≈4 at faint UV magnitudes. This may also have an effect on the difference in X_ LAE at z≈4 shown in Figure <ref>. Indeed, the difference for the 55 Å cut is more pronounced than that for the 25 Å cut. 
Although we do not discuss quantitatively biases of dropout selections at other redshifts, strong Lyα emission will cause similar effects, as discussed in the references. We note that the LBG sample in <cit.> consists of ≈70%  photometric-redshift selected objects and of ≈30% dropout selected objects based on dropout selection criteria for their medium-band filters. Since we measure X_ LAE for a photometric-redshift selected sample, this may result in the similarity of X_ LAE to ours at z≲5, though their sample is not complete in UV at M_1500≳-20 mag. 

Note that <cit.> mention yet another potential bias for LBG samples which are incomplete in UV, due to a correlation between the EW( Lyα) and M_1500. This bias leads to an underestimate of X_ LAE because large equivalent width objects are preferentially missed when faint-UV galaxies drop out of the sample. Our results are not affected by this bias, but it could affect the results of other work shown in Figure <ref>, and could have compensated the LBG selection bias we discussed above. In Figure <ref>, the bias from <cit.> has no impact because we are looking at X_ LAE as a function of UV magnitude. 






  §.§.§ Different observational methods


The Lyα emission in our sample is measured with an IFU (without including the Lyα halo) and is thus less affected by uncertainties due to slit-loss and aperture corrections. <cit.> measure the spatial offset between the Lyα emission and the UV continuum. They find a typical standard deviation for the offset which decreases towards higher redshifts (2.17^+0.19_-0.14 kpc (≈03) at z≈3.25 to 1.19^+1.29_-0.33 kpc (≈02) at z≈5.25). They argue that the evolution of the spatial offset contributes to the increasing trend of X_ LAE with z measured with slit spectroscopy with 1 slits such as in <cit.>. According to <cit.> Figure 7, the simulated cumulative distribution function (CDF) of slit-loss is similar from z≈3.5 to 5.5 but is shifted at z≈3–3.5 to a larger slit-loss for a 1 slit. At z≈3.5-4.5, the CDF reaches ≈90% at a slitloss of ≈10%. Moreover, their measured offsets are much larger than that for a lensed LAE at z≈1.8 <cit.> and typical values for LAEs at z≈3–6 <cit.>. Hence, this is probably not the dominant cause of the high X_ LAE values of <cit.> at z≈3.5–4.5. Note that the aperture diameters (convolved mask diameters) for our Lyα measurement with IFU data are typically larger than 1” (see Sect. <ref>). Our measurements are less affected by the spatial offset between the Lyα emission and the UV continuum. Meanwhile, <cit.> estimate slitlossess based on the spatial component of slit spectra. 




  §.§.§ Summary


We find indications that the dominant cause of the difference in X_ LAE measured in <cit.> and presented here is the LBG selection bias. Strong Lyα emission can enhance the probability to meet the LBG-selection criteria both in terms of the signal-to-noise ratio and color. The LBG selection bias has a strong effect on X_ LAE especially for faint UV magnitudes, where LBG samples are not complete. Possible discrepancies arising from different observational methods probably affect X_ LAE to a lesser extent. Thanks to the MUSE observations and to the HST photo-z sample, our X_ LAE are derived from the most homogeneous and complete sample to date.







 §.§ Comparison with the GALICS model

Using a homogeneous and complete UV sample and MUSE spectroscopic data, we have measured the LAE fraction for the first time at very faint magnitudes (M_1500≤-17.0 mag). While we confirm a weak increase of X_ LAE as a function of redshift at 3 ≲ z ≲ 5, we find that LAEs with EW( Lyα)≥45 Å  make up a relatively low fraction of the underlying rest-frame UV-detected galaxy population, ≈0–20%. This implies the existence of a duty cycle either for the star formation activity or for the escape and/or production of Lyα photons. Another possibility is that only a small fraction of all galaxies can evolve into LAEs or can be observed as LAEs in a limited range of inclinations <cit.>. Our results suggest no dependence of X_ LAE on M_1500 at z≈3. Keeping in mind that we want to assess the merits of the redshift evolution of X_ LAE to probe the IGM neutral fraction at z≳ 6, it is essential to understand these trends after reionization. To do so, we compare our results with predictions from the semi-analytic model of <cit.>. 



  §.§.§ Description of the model


<cit.> present an updated version of the GALICS hybrid model <cit.> which is designed to study the formation and evolution of galaxies in the high redshift Universe. GALICS relies on an N-body cosmological simulation to follow the hierarchical growth of dark matter structures and on semi-analytic prescriptions to describe the physics of the baryonic component. The box size of the simulation is 100 h^-1 cMpc on a side, and the dark-matter particle mass is ≈8.5×10^7  M_⊙ (with 1024^3 particles) <cit.>. In <cit.>, stars are formed according to a Kennicutt-Schmidt law when the galaxy’s gas surface density Σ_ gas is larger than a threshold value, Σ_ gas^ thresh <cit.>, and the intrinsic  emission from galaxies is computed assuming case B recombination as L_Lyα^ intr∝ 0.67 Q(H). Here, Q(H) is the production rate of hydrogen-ionising photons estimated from the stellar spectral energy distributions.
In order to predict the observed  properties of galaxies, <cit.> combine GALICS with the library of radiative transfer (RT) simulations of <cit.> which predict the escape fraction of  photons through galactic winds <cit.>. f_ esc depends on the wind parameters (the wind expansion velocity, velocity dispersion, dust opacity, neutral hydrogen column density) which are computed by GALICS. The  luminosity emerging from each individual galaxy is then given by L_ Lyα = L_ Lyα^ intr× f_ esc. 

The GALICS model was tuned to reproduce the UV and Lyα luminosity functions at z≈3–6 in <cit.>. This model was then shown to also reproduce accurately the observed stellar mass functions and star-formation-rate to stellar mass relations at z≈3–6 <cit.>. While their fiducial model can match these observational constraints at 3 ≲ z ≲ 6, it fails to reproduce the wide distribution of  EWs, in particular the high EW values (at EW( Lyα)≳50 Å). <cit.> discuss the possibility that this mismatch is linked to the lack of burstiness of star formation in GALICS given that the  EW is primarily set by the combination of (i) the production rate of  photons, dominated by short-lived stars, and (ii) the stellar UV emission which traces star formation over longer timescales <cit.>. In the fiducial model, most galaxies keep forming stars at a rather constant rate because the surface gas density threshold is almost always met. In an alternative model (labeled bursty SF), they increase this threshold by a factor 10, such that gas needs to accrete onto galaxies for longer periods before reaching the required surface density. This naturally gives rise to a star formation duty cycle and <cit.> show that their bursty SF model predicts EW distributions in much better agreement with observations than the fiducial model does. 

Following the procedure of <cit.>, we create mock surveys for both the fiducial and bursty SF models using the Mock Map Facility (MOMAF) tool <cit.>. In practice, for each model, we generate 100 lightcones that mimic the geometry and redshift range of the MUSE HUDF survey, i.e. a square field of ≈10 arcmin^2 and 2.8 ≲ z ≲ 6.7, and we compute X_ LAE from the mocks in the same bins of redshift and UV magnitude as for the observational measurements.  



  §.§.§ Measured LAE fraction vs. GALICS predictions


In Figure <ref>, we show our MUSE measurement of X_ LAE as a function of M_1500 at z≈3.3 for EW( Lyα)≥45 Å and EW( Lyα)≥65 Å. These EW cuts correspond to our most secure measurements. We also show in this Figure the predictions from the GALICS model. 
Both fiducial and bursty SF GALICS models (dashed and solid lines, respectively) show an increase of X_ LAE towards faint UV magnitudes with a slope which is in good agreement with the data (star symbols). As discussed in <cit.> and <cit.>, this trend may be the result of two factors. First, UV bright sources have intrinsically smaller EW( Lyα) due to less significant or less recent bursts of star formation. Second, these galaxies are often more massive with higher  and dust contents which can dramatically reduce the escape fraction of Lyα photons and therefore the observed EW. Hence, few bright UV galaxies display a strong  emission line. 

For a more detailed comparison, we see that the fiducial GALICS model (dashed lines) does not reproduce the observed X_ LAE. This model overestimates (underestimates) the LAE fraction at almost all UV magnitudes for EW( Lyα)≥45 Å (EW( Lyα)≥65 Å). This is a consequence of the too narrow EW distribution predicted by this model (see Sect. <ref>), which <cit.> attribute to overly smooth star formation histories. In the bursty SF model however (solid lines), galaxies have more diverse recent star formation histories which result in wider EW distributions, and consequently this model is able to reproduce our measured X_ LAE much better. 

In Figure <ref>, we compare the GALICS predictions of X_ LAE as a function of z with X_ LAE from the MUSE data at z≲5, i.e. where our observational measurements are most robust. For the same reasons as above, we find that the bursty SF model is much more successful at reproducing the observations than the fiducial model. This is particularly true for the lowest EW cut (i.e. 45 Å) where the agreement is quite good (solid orange curve). For EW( Lyα)≥65 Å however, we note that the bursty model slightly underpredicts the observed LAE fraction (solid magenta curve), especially at z ≳ 4.5. Additional ingredients could possibly be missing from this model that would help produce more galaxies with large EWs (in particular in the higher redshift bin) such as radiative transfer in asymmetric geometries, or Lyα production from other channels like collisions (gravitational cooling) or fluorescence <cit.>. Also, the assumed IMF or the metallicity evolution of model galaxies may not be realistic and lead to low EWs <cit.>.


Overall, these comparisons suggest that the measurements of X_ LAE by MUSE in the post-reionization epoch can be reasonably well interpreted with current models of high-z galaxies such as GALICS. We find that the observed trends between X_ LAE and redshift/UV magnitude are mainly shaped by the burstiness of star formation in GALICS. It is also caused by the variation of f_ esc with respect to the physical properties of the galaxies as discussed in <cit.>. In GALICS, these two aspects modulate the observed  EWs of galaxies and therefore the LAE fraction at z ≲ 5.














  §.§.§ Cosmic variance




The area of the MUSE HUDF survey is limited to 9.92 arcmin^2 which translates into comoving volumes of ≈ 1.5-2.5 × 10^4 cMpc^3 for the redshift ranges we are considering here. As explained in Sect. <ref>, we have accounted for several sources of uncertainty to compute the error on X_ LAE but so far we ignored cosmic variance (see Sect. <ref> for discussion). To assess the significance of this effect, we can estimate the cosmic variance from the GALICS mock lightcones, which are cut out from the ≈3×10^6 cMpc^3 simulation box <cit.>. We compute the 1σ standard deviation as the field-to-field variation, which includes the effects of both cosmic variance and the binomial proportion confidence interval. Note that our estimate of cosmic variance with GALICS only accounts for the clustering of galaxies and not for the possible contribution of large-scale variations in the IGM transmissivity due to a patchy reionization. 

In the middle and bottom panels of Figure <ref>, we compare the relative uncertainty of our MUSE X_ LAE estimates (circles) with the relative uncertainties due to the field-to-field variation of X_ LAE (solid lines with dots), which is estimated from the 100 mocks based on the bursty SF model. For a fair comparison, we match X_ LAE of GALICS to that of MUSE for EW( Lyα)≥45 Å (EW( Lyα)≥65 Å), by adopting slightly different EW( Lyα) cuts in the model catalogs, of 46 Å, 48 Å, and 46Å (53 Å, 52 Å, and 49 Å) at z≈3.3, 4.1, and 4.7,  respectively (see the top panel of Figure <ref>). Note that the total uncertainty for the MUSE X_ LAE is calculated by summing the statistical error (binomial proportional confidence interval) multiplied by completeness corrections in each flux bin for completeness correction (see Sect. <ref> for more details). 
For both EW cuts, the relative upper errors of our MUSE X_ LAE are much larger than those of the field-to-field variance for the bursty GALICS model. The relative lower errors of our MUSE X_ LAE are in the same level of those of the field-to-field variance at z≈4–5. The contribution of cosmic variance to the field-to-field variance in the GALICS mock is negligible, since statistical errors (crosses) are dominant. It suggests that the cosmic variance is a subdominant source of uncertainties in our measurement of X_ LAE. Therefore, we conclude that our MUSE results are not strongly affected by cosmic variance. We note that the field-to-field variance may be slightly underestimated in the mock catalogs because the fluctuations on scales larger than the simulated box (150 cMpc) are not sampled. The size of the simulated volume is however significantly larger than our MUSE survey volume (≈ 1.5-2.5 × 10^4 cMpc^3) and so this underestimate should be weak.







 §.§ The redshift evolution of the LAE fraction and implications for reionization



The main purpose of this paper is to measure the evolution of X_ LAE in the post-reionization epoch, at z≲6 (as shown in Figures <ref>, <ref>, and <ref>). Our results confirm the rise of X_ LAE with redshift found in the literature between z≈3 and 6 for -21.75≤ M_1500≤-17.75 mag, X_ LAE = 0.07^+0.06_-0.03z -0.22^+0.12_-0.24, in Figure <ref>. Meanwhile, the trend stopps at z∼5.6 for -20.25≤ M_1500≤-18.75 mag in Figure <ref>.  As discussed in Sect. <ref>, this evolution at z≈3-5 is not caused by the cosmic variance of the limited survey field of the MUSE HUDF Survey. Instead, it is probably caused by higher intrinsic EW( Lyα) and/or higher Lyα escape fractions at higher redshift in a given M_1500 range, due to less massive and less dusty galaxies at higher redshift <cit.>. It is also very important to understand the co-evolution of the Lyα and UV luminosity functions <cit.>. 

<cit.> obtain a relatively low X_ LAE (≈40%) at z≈6, which implies a less dramatic turn-over at z>6 than previously found <cit.>. If  we interpret our point at z≈ 5.5 as a statistical fluctuation 1σ below a true value ≈ 0.35 as found by <cit.>, we confirm this shallower increase of the LAE fraction towards z≈5 and z≈6. This could indicate the stop of the evolution of the Lyα escape fraction, possibly related to the plateau evolution of the star formation main sequence at z≈5–6 suggested by <cit.> and <cit.>, and implying a constant stellar mass at a given M_1500 <cit.>. 


Another possibility is that our low value of X_ LAE at z≈5.5 is genuine and indeed indicates a late transition in the ionisation state of the IGM. However, our data at z≤5 combined the work of <cit.> and <cit.>, which is not affected either by the LBG selection bias discussed above, suggest an earlier reionization, at z≈6-7. Our measurement at z≈5.5 may thus indicate a patchy reionization process. <cit.> measure the mean and scatter of the IGM Lyα opacity with the largest sample of quasars so far. They confirm the existence of tails towards high values in the Lyα opacity distributions, which may persist down to z≈5.2. They find a linear increase in the mean Lyα opacity from ≈1.8 at z≈5 to ≈3.8 at z≈6. These results also imply a late and/or patchy reionization scenario, in which reionization ends at z≈5.2–5.3 <cit.>. The Gunn-Peterson absorption trough in quasar spectra is only sensitive to a low Lyα opacity (very low  gas fraction) and the LAE fraction is therefore a complementary tool. In the near future, the James Web Space Telescope (JWST)/Near Infrared Spectrograph (NIRspec)  will enable us to observe Lyα emission at z≈5 to z≳10 homogeneously and help make significant progress. We can also use the WST/Near Infrared Imager and Slitless Spectrograph (NIRISS) for Lyα spectroscopy. Most importantly, one can measure Hα emission at z≈0 to z≈7 with JWST/NIRspec, and subsequently the line ratio of Lyα to Hα to disentangle between intrinsic evolution and escape fraction.

Another important point raised in this paper is that X_ LAE estimates are sensitive to EW( Lyα) selections (see Figures <ref> and <ref>). Although a general consensus seems to emerge from previous work, a quantitative interpretation of the evolution of X_ LAE with redshift requires more accurately constructed samples. In addition, the contribution of extended Lyα emission to the total Lyα budget is typically large and has a large scatter <cit.>. Methods for measuring Lyα flux have a large effect on EW( Lyα) and then X_ LAE. It means that accurate and homogeneous measurements of Lyα emission are required to use X_ LAE as a tracer of the  gas fraction of the IGM. In Figures <ref> and <ref>, the X_ LAE-z relation does not depend strongly on the rest-frame UV absolute magnitude. This suggests that at z<6, combining UV-bright and faint samples can give us better statistics for X_ LAE measurements. In addition, as discussed in Sect. <ref>, a firm definition of parent samples avoiding a selection bias is also required to assess the evolution of X_ LAE. The uncertainties of X_ LAE have to be calculated with BPCI (Bernoulli trials, see Sect. <ref>). Moreover, <cit.> warn of the interpretation of the evolution of X_ LAE with the same UV magnitude, since galaxies with the same UV magnitude have very different stellar and halo masses at different redshifts <cit.>. Because of such effects, sophisticated models of galaxy formation are needed to robustly interpret variations of X_ LAE with cosmic time. 
We propose a method using a UV complete sample including faint galaxies, based on a photo-z selection and an absolute magnitude cut, together with Lyα measurements by an IFU with a high sensitivity and a wide-wavelength coverage in a large field-of-view like VLT/MUSE and VLT/BlueMUSE <cit.> at z≈2–6.6. 








 
 














§ SUMMARY AND CONCLUSIONS

We have investigated the LAE fraction at z≈3–6 using the second data release of the  MUSE Hubble Ultra Deep Field Survey and the HST catalog of the UVUDF. Thanks to the unprecedented depth of the MUSE and HST data for Lyα and UV, respectively, we have studied the LAE fraction for galaxies as faint as M_1500=-17.0 mag at z≈3 for the first time with a UV-complete sample. We have also derived the LAE fraction as a function of redshift homogeneously from z=3 to 6, down to M_1500=-17.75 mag. Our results are summarized as follows: 



  * We derived the redshift evolution of X_ LAE for a number of EW and UV magnitude selections, including the first estimate down to -17.75 mag. These results are summarised in Table <ref>. For all selections, we find low values of X_ LAE ≈0.04-0.3, and a weak rise of X_ LAE with z, qualitatively consistent with the trend reported for brighter samples in the literature.



  * We compared our MUSE results with those in the literature for M_1500∈ [-20.25;-18.75]. At z≲5, our values of X_ LAE are consistent with those in <cit.> and <cit.> within 1σ error bars for EW( Lyα)≥25 Å (see left panel of Figure <ref>). Our X_ LAE at z≈5.6 is lower than those in the literature, which may be caused by statistical errors, or a late and/or patchy reionization process.




  * We measured the dependence of X_ LAE on M_1500 at z=2.9–3.7 for EW( Lyα)≥25 Å, 45 Å, 65 Å, and 85 Å (see Figure <ref>). This is the first time this has been measured down to M_1500=-17.0 mag (for the largest EWs of our sample), and for a volume-limited sample. We found no clear dependence of X_ LAE on M_1500, in contrast to previous reports. 


  * We compared the dependence of X_ LAE on M_1500 for EW( Lyα)≥50 Å at z≈3-4 derived from MUSE with results from the literature (Figure <ref>). Again we found no dependence of X_ LAE on M_1500. Our slopes of 0.01^+0.02_-0.03 and 0.01^+0.06_-0.10 at z≈2.9–3.7 and at ≈3.7–4.4, respectively, are shallower than that in <cit.>, 0.13^+0.03_-0.03 at ≈3.5–4.5. We also found lower values of our X_ LAE at a faint UV magnitude of M_1500≳-19 mag. 


  * The dominant causes of the difference of X_ LAE in our work and previous studies appear to be LBG selection biases in those studies. We showed how these can lead to an over-estimate by a factor ≈1.5 of X_ LAE at z≈4 for galaxies with M_1500=-19 mag and EW( Lyα)≥50 Å.


  * We compared our MUSE results with predictions from a cosmological semi-analytic galaxy evolution model <cit.>. When GALICS uses a bursty star formation model, it can reproduce our measurement of X_ LAE as a function of M_1500 at z≈3. The fiducial GALICS model however cannot. The bursty model can also reproduce X_ LAE as a function of z at z≲4. We assessed cosmic variance for our MUSE results using the bursty SF model and found that it does not have a significant effect on our results.


  * Overall, we found that X_ LAE is lower than ≈30%. This implies a low duty cycle of LAEs, suggesting bursty star formation or strong time variations in the production of Lyα photons and/or in their escape fraction. 





Despite the difficulties of the method, the dominant source of uncertainties in our work is the Poisson noise due to the small number of objects in our samples. This is encouraging and suggests that future deep surveys with e.g., MUSE and JWST will enable us to produce accurate measurements of X_ LAE with secure samples and to extend our understanding of the evolution of X_ LAE at all redshifts, after and during the epoch of reionization. 



We thank the anonymous referee for constructive comments and suggestions. We would like to express our gratitude to Stephane De Barros and Pablo Arrabal Haro for kindly providing their data plotted in Figures <ref>, <ref>,  and <ref>. We are grateful to Pascal Oesch, Kazuhiro Shimasaku, Masami Ouchi, Rieko Momose, Daniel Schaerer, Hidenobu Yajima, Taku Okamura, Makoto Ando, and Hinako Goto for giving insightful comments and suggestions. This work is based on observations taken by VLT, which is operated by European Southern Observatory. This research made use of Astropy[<http://www.astropy.org>], which is a community-developed core Python package for Astronomy <cit.>, , , and  <cit.>. H.K. acknowledges support from Japan Society for the Promotion of Science (JSPS) through the JSPS Research Fellowship for Young Scientists and Overseas Challenge Program for Young Researchers. This work was supported by the project FOGHAR (Agence Nationale de la Recherche, ANR-13-BS05-0010-02). JB acknowledges support from the ORAGE project from the Agence Nationale de la Recherche under grant ANR-14-CE33-0016-03. JR acknowledges support from the ERC starting grant 336736-CALENDS. T. H. acknowledges supports by the Grant-inAid for Scientic Research 19J01620. 








§ UNCERTAINTIES OF Z_ P



 §.§ Impact of lacking IRAC data on z_ p estimation in <cit.>

It is well known that LBG samples and photo-z samples at z≈3–7 can be contaminated by lower-z galaxies with a 3646 Å Balmer or 4000 Å break at z≈0–1, especially at faint magnitude. Spitzer/IRAC data can provide rest-frame optical data which are useful to break the degeneracy of z_ p <cit.>. In this work, we use the catalog from R15, where z_ p are derived using HST data alone. The advantage of this choice is discussed in <cit.>, where they show that IRAC data can in fact worsen photo-z performance for faint galaxies in the MUSE HUDF sample (see their Appendix A). It might be a reflection of the difficulty of providing reliable IRAC photometry for sources as faint as most of our sample.

As discussed in Sect. <ref>, the fraction of low-z contaminants among R15 galaxies is found to be low within MUSE samples <cit.>. To avoid contaminants due to poor photo-z estimations, we apply an S/N>2 cut for our sample and then apply a stricter cut on M_1500 (see Sect. <ref> and Figure <ref> for more details).




 §.§ Effects of z_ p errors on redshift binning


To check the effect of the error on z_ p on redshift binning, we compare the median of upper and lower 95% errors of z_ p to the half-width of redshift bins for -20.25≤ M_1500≤-18.75 mag and -18.75≤ M_1500≤-17.75 mag shown in Figures <ref> and <ref>. The upper (lower) 95% errors are calculated from a difference between 95% upper (lower) limit of photo-z and photo-z where likelihood is maximized in BPZ. The results are summarized in Table <ref>. For -20.25≤ M_1500≤-18.75 mag, the median values of the 95% z_ p errors are smaller than the width of redshift bins at z≈3.3 to 5.6. For -18.75≤ M_1500≤-17.75 mag, the errors are larger than those for the brighter M_1500, but the median of 95% errors are still smaller than or comparable to the width of redshift bins at the all over redshift range. Therefore, the 95% errors of z_ p are typically smaller than the width of redshift bins in Figures <ref> and <ref>. Note that we include widths of the redshift bins in the linear relation fitting (see Sect. <ref>).

Next, we check the fraction of possible low-z contaminants, f_ large error. Although the probability distribution functions (PDFs) of z_ p in R15 catalog are not published, galaxies with a bimodal PDF of z_ p show a large lower 95% error, which reaches at z≈0–1 for a sample at z≈3–7. We calculate f_ large error from the wavelengths of breaks, z_ p, and its 95% errors and summarize the results in Table <ref>. 
For -20.25≤ M_1500≤-18.75 mag, f_ large error is 0.02 to 0.09, implying that our X_ LAE in Figure <ref> is not affected significantly by contaminants. Meanwhile, for -18.75≤ M_1500≤-17.75 mag, f_ large error is 0.05 at z≈3.3 so that our X_ LAE is not suppressed significantly by contaminants in Figure <ref>. At z≈4.1 to 5.6, f_ large error is relatively high, 0.13 to 0.35. However, these are conservative estimations of the upper limit of low-z contaminated fraction, since not all of the galaxies with a large photo-z error locate at z≈0–1. In fact, our sample shows the maximum likelihood at z≈3–6. Our X_ LAE at z≳4.1 and 5.6 in Figure <ref> should not be affected significantly by low-z contaminants significantly. 


 







§ THE TEMPLATE SPECTRA USED IN 
 
Figure <ref> shows the template spectra of LAEs in  used in this work. The continuum are subtracted in the templates and the bright Lyα lines are clipped to lie between -30 and +30 times the mean absolute deviation in the similar manner to those used in the original  <cit.> and  <cit.> as shown in the right panels in Figure <ref>. Cross-correlation functions indicate locations of lines in the fitting range and are not affected significantly by the line shape of templates in general <cit.>. In fact, the completeness of  does not depend on the FWHM of fake lines in completeness simulations as we mention in Sect. <ref>. Although our templates do not cover various types of Lyα lines, it does not have an effect on detection of Lyα emission.  







§ AN EXAMPLE OF CONTAMINATION OF LYΑ DETECTION

Figure <ref> shows an example of contamination of Lyα emission from a neighboring object. Panels (a) - (d) show sub-panels in  for a UV-selected source <cit.>. The HST cutout around a UV-galaxy with <cit.> ID =628 is shown in panel (a). The 1D spectrum shown in panel (d) clearly shows a strong Lyα emission line with the highest confidence level by . The spectrum is extracted from the object mask (panel (b)), and clearly contaminated by diffuse Lyα emission from a neighboring object (MUSE ID =1185 in the DR1 catalog) as shown in the narrow band in panel (c). We can remove these contaminated objects from our sample of Lyα emitter candidates in visual inspection.  






§ THE NUMBER OF FLUX BINS USED TO CORRECT INCOMPLETENESS OF LYΑ DETECTION


We also examine the effect of binning to correct incompleteness of the number of LAEs. Because the number of objects we have is small, somewhat arbitrary binning may affect X_ LAE. For a large number of bins, we often get flux-bins with no object at all, and these bins increase the error bars. For a small number of bins, we introduce a large error on the completeness correction as described above. To test the effect of binning, we vary the number of bins (N_ bin) we use, from 2 to 6, and see how the median values and error bars change. In Figure <ref>, we show the median values and error bars of X_ LAE for N_ bin=2–6 for the plot of the evolution of X_ LAE (Figure<ref>). The uncertainties of the completeness correction are 20%, 25%, 32%, 44%, and 73% for N_ bin=6, 5, 4, 3, and 2 (see Sect. <ref>).  We find that 4 to 6 bins are a sweet spot where the error bars are small and appear converged and adopt N_ bin=4 in Sect. <ref>. 









§ ERROR PROPAGATION OF COMPLETENESS CORRECTION VALUES IN A BINOMIAL PROPORTION CONFIDENCE INTERVAL




We test the applicability of a binomial proportion confidence interval (BPCI) for the case of a completeness correction to calculate error bars of the LAE fraction. First, we examine how to propagate completeness correction values in the error calculation with BPCI numerically. To do mock observations, we generate N_ LAE^ true(z_ s, M_1500,  EW) randomly with 100000 trials using the python module  for each N_ 1500(z_ p,  M_1500) and the true LAE fraction (X_ LAE^ true). Then we generate N_ LAE^ det(z_ s, M_1500,  EW) randomly for each N_ LAE^ true(z_ s, M_1500,  EW) again using  for a given completeness correction value as a probability of detection. We can obtain the probability distribution of X_ LAE with a given completeness correction value, N_ 1500(z_ p,  M_1500), and X_ LAE^ true. We compare the 1σ upper and lower uncertainties from the probability distribution function of the mock X_ LAE with those derived from two methods of BPCI () with input parameters of the number of LAEs and N_ 1500(z_ p,  M_1500) if the incompleteness of the number of LAEs are corrected (N_ LAE^ true(z_ s, M_1500,  EW)) or not (N_ LAE^ det(z_ s, M_1500,  EW)). We confirm that it is better to input N_ LAE^ det(z_ s, M_1500,  EW) and to multiply the obtained uncertainties by the correction value (see also Sect. <ref>). 

Second, we check the accuracy of the method above in the plane of N_ 1500(z_ p,  M_1500) and X_ LAE^ true for a given completeness correction value. We calculate the 1σ upper and lower limits of X_ LAE with this method. Again we generate 100000 mock values of N_ LAE^ det(z_ s, M_1500,  EW) and then X_ LAE numerically. The fractions of the mock X_ LAE within the range of the 1σ upper and lower limits among all the experiments are calculated. If the method is accurate, this fraction would be ≈68%. We check the fraction of experiments for N_ 1500(z_ p,  M_1500) from 0 to 250 with a step of 1 and for X_ LAE^ true from 0.1 to 0.5 with a step of 0.02. In Figure <ref>, the example with low completeness correction values, 0.1 and 0.5, is shown. In the panel (a) for completeness correction =0.1, most of the plane is colored with yellow green corresponding to a fraction close to 0.68, while it is colored with blue or red corresponding to overestimation and underestimation of the errors, respectively, for cases with poor statistics (i.e., low N_ 1500(z_ p,  M_1500) and X_ LAE^ true). In panel (b) for completeness correction =0.5, the method is found to produce errors accurately except for the very poor statistics cases (at the upper left region in the panel). We can estimate the error more accurately with this method for a higher completeness case. Even with a low completeness and smallest numbers, the uncertainties are overestimated. Therefore, we adopt this method conservatively.  






§ THE BEST FIT LINEAR RELATIONS OF X_ LAE AS A FUNCTION OF Z AND M_1500

 We show the best fit linear relations of X_ LAE as a function of z and M_1500 in Figure <ref>. The equations of relations are shown in Sects. <ref> and <ref>. 





aa
126
natexlab#1#1[Arrabal Haro et al.(2018)Arrabal Haro, Rodríguez Espinosa,
  Muñoz-Tuñón, Pérez-González, Dannerbauer,
  Bongiovanni, Barro, Cava, Lumbreras-Calle, Hernán-Caballero,
  Eliche-Moral, Domínguez Sánchez, Conselice, Tresse,
  Alcalde Pampliega, Balcells, Daddi, & Rodighiero]ArrabalHaro2018
Arrabal Haro, P., Rodríguez Espinosa, J. M.,
  Muñoz-Tuñón, C., et al. 2018, Monthly Notices of the
  Royal Astronomical Society, 478, 3740

[Bacon et al.(2010)Bacon, Accardo, Adjali, Anwand, Bauer, Biswas,
  Blaizot, Boudon, Brau-Nogue, Brinchmann, Caillier, Capoani, Carollo, Contini,
  Couderc, Daguisé, Deiries, Delabre, Dreizler, Dubois, Dupieux, Dupuy,
  Emsellem, Fechner, Fleischmann, François, Gallou, Gharsa, Glindemann,
  Gojak, Guiderdoni, Hansali, Hahn, Jarno, Kelz, Koehler, Kosmalski, Laurent,
  Le Floch, Lilly, Lizon, Loupias, Manescau, Monstein, Nicklas, Olaya, Pares,
  Pasquini, Pécontal-Rousset, Pelló, Petit, Popow, Reiss,
  Remillieux, Renault, Roth, Rupprecht, Serre, Schaye, Soucail, Steinmetz,
  Streicher, Stuik, Valentin, Vernet, Weilbacher, Wisotzki, &
  Yerle]Bacon2010
Bacon, R., Accardo, M., Adjali, L., et al. 2010, in Proceedings of the SPIE,
  Volume 7735, ed. I. S. McLean, S. K. Ramsay, & H. Takami, 773508

[Bacon et al.(2015)Bacon, Brinchmann, Richard, Contini, Drake,
  Franx, Tacchella, Vernet, Wisotzki, Blaizot, Bouché, Bouwens,
  Cantalupo, Carollo, Carton, Caruana, Clément, Dreizler, Epinat,
  Guiderdoni, Herenz, Husser, Kamann, Kerutt, Kollatschny, Krajnovic, Lilly,
  Martinsson, Michel-Dansac, Patricio, Schaye, Shirazi, Soto, Soucail,
  Steinmetz, Urrutia, Weilbacher, & de Zeeuw]Bacon2015
Bacon, R., Brinchmann, J., Richard, J., et al. 2015, Astronomy &
  Astrophysics, 575, A75

[Bacon et al.(2017)Bacon, Conseil, Mary, Brinchmann, Shepherd,
  Akhlaghi, Weilbacher, Piqueras, Wisotzki, Lagattuta, Epinat, Guerou, Inami,
  Cantalupo, Courbot, Contini, Richard, Maseda, Bouwens, Bouché,
  Kollatschny, Schaye, & Pello]Bacon2017
Bacon, R., Conseil, S., Mary, D., et al. 2017, Astronomy & Astrophysics,
  608

[Baldry et al.(2014)Baldry, Alpaslan, Bauer, Bland-Hawthorn, Brough,
  Cluver, Croom, Davies, Driver, Gunawardhana, Holwerda, Hopkins, Kelvin,
  Liske, Lopez-Sanchez, Loveday, Norberg, Peacock, Robotham, &
  Taylor]Baldry2014
Baldry, I. K., Alpaslan, M., Bauer, A. E., et al. 2014, Monthly Notices of
  the Royal Astronomical Society, 441, 2440

[Behroozi et al.(2013)Behroozi, Wechsler, & Conroy]Behroozi2013
Behroozi, P. S., Wechsler, R. H., & Conroy, C. 2013, The Astrophysical
  Journal, 770, 57

[Benitez(2000)]Benitez2000
Benitez, N. 2000, The Astrophysical Journal, 536, 571

[Benitez et al.(2004)Benitez, Ford, Bouwens, Menanteau, Blakeslee,
  Gronwall, Illingworth, Meurer, Broadhurst, Clampin, Franx, Hartig, Magee,
  Sirianni, Ardila, Bartko, Brown, Burrows, Cheng, Cross, Feldman, Golimowski,
  Infante, Kimble, Krist, Lesser, Levay, Martel, Miley, Postman, Rosati,
  Sparks, Tran, Tsvetanov, White, & Zheng]Benitez2004
Benitez, N., Ford, H., Bouwens, R., et al. 2004, The Astrophysical Journal
  Supplement Series, 150, 1

[Blaizot et al.(2005)Blaizot, Wadadekar, Guiderdoni, Colombi,
  Bertin, Bouchet, Devriendt, & Hatton]Blaizot2005
Blaizot, J., Wadadekar, Y., Guiderdoni, B., et al. 2005, Monthly Notices of
  the Royal Astronomical Society, 360, 159

[Bosman et al.(2018)Bosman, Fan, Jiang, Reed, Matsuoka, Becker, &
  Haehnelt]Bosman2018
Bosman, S. E., Fan, X., Jiang, L., et al. 2018, Monthly Notices of the Royal
  Astronomical Society, 479, 1055

[Bouwens et al.(2016)Bouwens, Aravena, Decarli, Walter, da Cunha,
  Labbé, E. Bauer, Bertoldi, Carilli, Chapman, Daddi, Hodge, J. Ivison,
  Karim, Le Fevre, Magnelli, Ota, Riechers, R. Smail, van der Werf, Weiss, Cox,
  Elbaz, Gonzalez-Lopez, Infante, Oesch, Wagg, & Wilkins]Bouwens2016a
Bouwens, R. J., Aravena, M., Decarli, R., et al. 2016, The Astrophysical
  Journal, 833, 72

[Bouwens et al.(2007)Bouwens, Illingworth, Franx, &
  Ford]Bouwens2007
Bouwens, R. J., Illingworth, G. D., Franx, M., & Ford, H. 2007, The
  Astrophysical Journal, 670, 928

[Bouwens et al.(2015a)Bouwens, Illingworth, Oesch,
  Caruana, Holwerda, Smit, & Wilkins]Bouwens2015a
Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015a,
  The Astrophysical Journal, 811, 140

[Bouwens et al.(2015b)Bouwens, Illingworth, Oesch,
  Trenti, Labbé, Bradley, Carollo, Van Dokkum, Gonzalez, Holwerda, Franx,
  Spitler, Smit, & Magee]Bouwens2015b
Bouwens, R. J., Illingworth, G. D., Oesch, P. A., et al. 2015b,
  Astrophysical Journal, 803, 1

[Bradač et al.(2019)Bradač, Huang, Fontana, Castellano,
  Merlin, Amorín, Hoag, Strait, Santini, Ryan, Casertano, Lemaux, Lubin,
  Schmidt, Schrabback, Treu, von der Linden, Mason, & Wang]Bradac2019
Bradač, M., Huang, K.-H., Fontana, A., et al. 2019, Monthly Notices of
  the Royal Astronomical Society, 489, 99

[Brammer et al.(2008)Brammer, van Dokkum, & Coppi]Brammer2008
Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, The Astrophysical
  Journal, 686, 1503

[Brinchmann et al.(2017)Brinchmann, Inami, Bacon, Contini, Maseda,
  Chevallard, Bouché, Boogaard, Carollo, Charlot, Kollatschny, Marino,
  Pello, Richard, Schaye, Verhamme, & Wisotzki]Brinchmann2017
Brinchmann, J., Inami, H., Bacon, R., et al. 2017, Astronomy &
  Astrophysics, 608, A3

[Cameron(2011)]Cameron2011
Cameron, E. 2011, Publications of the Astronomical Society of Australia, 28,
  128

[Caruana et al.(2014)Caruana, Bunker, Wilkins, Stanway, Lorenzoni,
  Jarvis, & Ebert]Caruana2014
Caruana, J., Bunker, A. J., Wilkins, S. M., et al. 2014, Monthly Notices of
  the Royal Astronomical Society, 443, 2831

[Caruana et al.(2018)Caruana, Wisotzki, Herenz, Kerutt, Urrutia,
  Schmidt, Bouwens, Brinchmann, Cantalupo, Carollo, Diener, Drake, Garel,
  Marino, Richard, Saust, Schaye, & Verhamme]Caruana2018
Caruana, J., Wisotzki, L., Herenz, E. C., et al. 2018, Monthly Notices of the
  Royal Astronomical Society, 473, 30

[Cassata et al.(2015)Cassata, Tasca, Le Fèvre, Lemaux,
  Garilli, Le Brun, Maccagni, Pentericci, Thomas, Vanzella, Zamorani, Zucca,
  Amorin, Bardelli, Capak, Cassarà, Castellano, Cimatti, Cuby, Cucciati,
  de la Torre, Durkalec, Fontana, Giavalisco, Grazian, Hathi, Ilbert, Moreau,
  Paltani, Ribeiro, Salvato, Schaerer, Scodeggio, Sommariva, Talia, Taniguchi,
  Tresse, Vergani, Wang, Charlot, Contini, Fotopoulou, Koekemoer,
  López-Sanjuan, Mellier, & Scoville]Cassata2015
Cassata, P., Tasca, L. A. M., Le Fèvre, O., et al. 2015, Astronomy &
  Astrophysics, 573, A24

[Charlot & Fall(1993)]Charlot1993
Charlot, S. & Fall, S. M. 1993, Astrophysical Journal v.415, 415, 580

[Coe et al.(2006)Coe, Benítez, Sánchez, Jee, Bouwens, &
  Ford]Coe2006
Coe, D., Benítez, N., Sánchez, S. F., et al. 2006, The
  Astronomical Journal, 132, 926

[Cooke et al.(2014)Cooke, Ryan-Weber, Garel, &
  Díaz]Cooke2014
Cooke, J., Ryan-Weber, E. V., Garel, T., & Díaz, C. G. 2014, Monthly
  Notices of the Royal Astronomical Society, 441, 837

[Curtis-Lake et al.(2012)Curtis-Lake, McLure, Pearce, Dunlop,
  Cirasuolo, Stark, Almaini, Bradshaw, Chuter, Foucaud, &
  Hartley]Curtis-Lake2012
Curtis-Lake, E., McLure, R. J., Pearce, H. J., et al. 2012, Monthly Notices
  of the Royal Astronomical Society, 422, 1425

[De Barros et al.(2017)De Barros, Pentericci, Vanzella, Castellano,
  Fontana, Grazian, Conselice, Yan, Koekemoer, Cristiani, Dickinson,
  Finkelstein, & Maiolino]DeBarros2017
De Barros, S., Pentericci, L., Vanzella, E., et al. 2017, Astronomy &
  Astrophysics, 608, A123

[Dijkstra(2014)]Dijkstra2014
Dijkstra, M. 2014, Publications of the Astronomical Society of Australia, 31

[Dijkstra(2017)]Dijkstra2017
Dijkstra, M. 2017, Lecture notes for the 46th Saas-Fee winterschool, 13

[Dijkstra et al.(2011)Dijkstra, Mesinger, & Wyithe]Dijkstra2011
Dijkstra, M., Mesinger, A., & Wyithe, J. S. B. 2011, Monthly Notices of the
  Royal Astronomical Society, 414, 2139

[Drake et al.(2017a)Drake, Garel, Wisotzki, Leclercq,
  Hashimoto, Richard, Bacon, Blaizot, Caruana, Conseil, Contini, Guiderdoni,
  Herenz, Inami, Lewis, Mahler, Marino, Pello, Schaye, Verhamme, Ventou, &
  Weilbacher]Drake2017b
Drake, A. B., Garel, T., Wisotzki, L., et al. 2017a, Astronomy
  & Astrophysics, 608, A6

[Drake et al.(2017b)Drake, Guiderdoni, Blaizot,
  Wisotzki, Herenz, Garel, Richard, Bacon, Bina, Cantalupo, Contini, den Brok,
  Hashimoto, Marino, Pelló, Schaye, & Schmidt]Drake2017a
Drake, A. B., Guiderdoni, B., Blaizot, J., et al. 2017b, Monthly
  Notices of the Royal Astronomical Society, 471, 267

[Dunlop(2013)]Dunlop2013
Dunlop, J. S. 2013, in Astrophysics and Space Science Library (Springer-Verlag
  Berlin Heidelberg), 223–292

[Erb et al.(2019)Erb, Berg, Auger, Kaplan, Brammer, &
  Pettini]Erb2019arXiv
Erb, D. K., Berg, D. A., Auger, M. W., et al. 2019, arXiv: 1907.11733

[Faisst et al.(2014)Faisst, Capak, Carollo, Scarlata, &
  Scoville]Faisst2014
Faisst, A. L., Capak, P., Carollo, C. M., Scarlata, C., & Scoville, N. 2014,
  The Astrophysical Journal, 788, 87

[Fan et al.(2006)Fan, Strauss, Becker, White, Gunn, Knapp, Richards,
  Schneider, Brinkmann, & Fukugita]Fan2006
Fan, X., Strauss, M. A., Becker, R. H., et al. 2006, The Astronomical
  Journal, 132, 117

[Finkelstein et al.(2015)Finkelstein, Ryan, Papovich, Dickinson,
  Song, Somerville, Ferguson, Salmon, Giavalisco, Koekemoer, Ashby, Behroozi,
  Castellano, Dunlop, Faber, Fazio, Fontana, Grogin, Hathi, Jaacks, Kocevski,
  Livermore, McLure, Merlin, Mobasher, Newman, Rafelski, Tilvi, &
  Willner]Finkelstein2015a
Finkelstein, S. L., Ryan, R. E., Papovich, C., et al. 2015, The Astrophysical
  Journal, 810, 71

[Fontana et al.(2010)Fontana, Vanzella, Pentericci, Castellano,
  Giavalisco, Grazian, Boutsia, Cristiani, Dickinson, Giallongo, Maiolino,
  Moorwood, & Santini]Fontana2010
Fontana, A., Vanzella, E., Pentericci, L., et al. 2010, The Astrophysical
  Journal, 725, L205

[Forero-Romero et al.(2012)Forero-Romero, Yepes, Gottlöber, &
  Prada]Forero-Romero2012
Forero-Romero, J. E., Yepes, G., Gottlöber, S., & Prada, F. 2012,
  Monthly Notices of the Royal Astronomical Society, 419, 952

[Garel et al.(2015)Garel, Blaizot, Guiderdoni, Michel-Dansac, Hayes,
  & Verhamme]Garel2015
Garel, T., Blaizot, J., Guiderdoni, B., et al. 2015, Monthly Notices of the
  Royal Astronomical Society, 450, 1279

[Garel et al.(2012)Garel, Blaizot, Guiderdoni, Schaerer, Verhamme,
  & Hayes]Garel2012
Garel, T., Blaizot, J., Guiderdoni, B., et al. 2012, Monthly Notices of the
  Royal Astronomical Society, 422, 310

[Garel et al.(2016)Garel, Guiderdoni, & Blaizot]Garel2016
Garel, T., Guiderdoni, B., & Blaizot, J. 2016, Monthly Notices of the Royal
  Astronomical Society, 455, 3436

[Gunn & Peterson(1965)]Gunn1965
Gunn, J. E. & Peterson, B. A. 1965, The Astrophysical Journal, 142, 1633

[Guo et al.(2013)Guo, Ferguson, Giavalisco, Barro, Willner, Ashby,
  Dahlen, Donley, Faber, Fontana, Galametz, Grazian, Huang, Kocevski,
  Koekemoer, Koo, McGrath, Peth, Salvato, Wuyts, Castellano, Cooray, Dickinson,
  Dunlop, Fazio, Gardner, Gawiser, Grogin, Hathi, Hsu, Lee, Lucas, Mobasher,
  Nandra, Newman, & van der Wel]Guo2013
Guo, Y., Ferguson, H. C., Giavalisco, M., et al. 2013, The Astrophysical
  Journal Supplement Series, 207, 24

[Hashimoto et al.(2017a)Hashimoto, Garel, Guiderdoni,
  Drake, Bacon, Blaizot, Richard, Leclercq, Inami, Verhamme, Bouwens,
  Brinchmann, Cantalupo, Carollo, Caruana, Herenz, Kerutt, Marino, Mitchell, &
  Schaye]Hashimoto2017b
Hashimoto, T., Garel, T., Guiderdoni, B., et al. 2017a,
  Astronomy & Astrophysics, 608, A10

[Hashimoto et al.(2017b)Hashimoto, Ouchi, Shimasaku,
  Schaerer, Nakajima, Shibuya, Ono, Rauch, & Goto]Hashimoto2017a
Hashimoto, T., Ouchi, M., Shimasaku, K., et al. 2017b, MNRAS,
  465, 1543

[Hashimoto et al.(2015)Hashimoto, Verhamme, Ouchi, Shimasaku,
  Schaerer, Nakajima, Shibuya, Rauch, Ono, & Goto]Hashimoto2015
Hashimoto, T., Verhamme, A., Ouchi, M., et al. 2015, The Astrophysical
  Journal, 812, 157

[Hatton et al.(2003)Hatton, Devriendt, Ninin, Bouchet, Guiderdoni,
  & Vibert]Hatton2003
Hatton, S., Devriendt, J. E. G., Ninin, S., et al. 2003, Monthly Notices of
  the Royal Astronomical Society, 343, 75

[Herenz & Wisotzki(2017)]Herenz2017b
Herenz, E. C. & Wisotzki, L. 2017, Astronomy & Astrophysics, 602, A111

[Herenz et al.(2019)Herenz, Wisotzki, Saust, Kerutt, Urrutia,
  Diener, Schmidt, Marino, de la Vieuville, Boogaard, Schaye, Guiderdoni,
  Richard, & Bacon]Herenz2019
Herenz, E. C., Wisotzki, L., Saust, R., et al. 2019, Astronomy &
  Astrophysics, 621, A107

[Hinton et al.(2016)Hinton, Davis, Lidman, Glazebrook, &
  Lewis]Hinton2016
Hinton, S. R., Davis, T. M., Lidman, C., Glazebrook, K., & Lewis, G. F. 2016,
  Astronomy and Computing, 15, 61

[Hoag et al.(2019a)Hoag, Bradač, Huang, Mason,
  Treu, Schmidt, Trenti, Strait, Lemaux, Finney, & Paddock]Hoag2019aApJ
Hoag, A., Bradač, M., Huang, K., et al. 2019a, The
  Astrophysical Journal, 878, 12

[Hoag et al.(2019b)Hoag, Treu, Pentericci, Amorin,
  Bolzonella, Bradač, Castellano, Cullen, Fynbo, Garilli, Guaita, Hathi,
  Henry, Jones, Mason, McLeod, McLure, Morishita, Pozzetti, Schaerer, Schmidt,
  Talia, & Thomas]Hoag2019bMNRAS
Hoag, A., Treu, T., Pentericci, L., et al. 2019b, Monthly
  Notices of the Royal Astronomical Society, 488, 706

[Hu et al.(2010)Hu, Cowie, Barger, Capak, Kakazu, &
  Trouille]Hu2010
Hu, E. M., Cowie, L. L., Barger, A. J., et al. 2010, The Astrophysical
  Journal, 725, 394

[Hunter(2007)]Hunter2007
Hunter, J. D. 2007, Computing in Science & Engineering, 9, 90

[Inami et al.(2017)Inami, Bacon, Brinchmann, Richard, Contini,
  Conseil, Hamer, Akhlaghi, Bouché, Clément, Desprez, Drake,
  Hashimoto, Leclercq, Maseda, Michel-Dansac, Paalvast, Tresse, Ventou,
  Kollatschny, Boogaard, Finley, Marino, Schaye, & Wisotzki]Inami2017
Inami, H., Bacon, R., Brinchmann, J., et al. 2017, Astronomy &
  Astrophysics, 608, A2

[Inoue et al.(2018)Inoue, Hasegawa, Ishiyama, Yajima, Shimizu,
  Umemura, Konno, Harikane, Shibuya, Ouchi, Shimasaku, Ono, Kusakabe, Higuchi,
  & Lee]Inoue2018
Inoue, A. K., Hasegawa, K., Ishiyama, T., et al. 2018, Publications of the
  Astronomical Society of Japan, 70, 1

[Itoh et al.(2018)Itoh, Ouchi, Zhang, Inoue, Mawatari, Shibuya,
  Harikane, Ono, Kusakabe, Shimasaku, Fujimoto, Iwata, Kajisawa, Kashikawa,
  Kawanomoto, Komiyama, Lee, Nagao, & Taniguchi]Itoh2018
Itoh, R., Ouchi, M., Zhang, H., et al. 2018, The Astrophysical Journal, 867,
  46

[Jensen et al.(2013)Jensen, Laursen, Mellema, Iliev, Sommer-Larsen,
  & Shapiro]Jensen2013
Jensen, H., Laursen, P., Mellema, G., et al. 2013, Monthly Notices of the
  Royal Astronomical Society, 428, 1366

[Kakiichi et al.(2016)Kakiichi, Dijkstra, Ciardi, &
  Graziani]Kakiichi2016
Kakiichi, K., Dijkstra, M., Ciardi, B., & Graziani, L. 2016, Monthly Notices
  of the Royal Astronomical Society, 463, 4019

[Kashikawa et al.(2006)Kashikawa, Shimasaku, Malkan, Doi, Matsuda,
  Ouchi, Taniguchi, Ly, Nagao, Iye, Motohara, Murayama, Murozono, Nariai, Ohta,
  Okamura, Sasaki, Shioya, & Umemura]Kashikawa2006
Kashikawa, N., Shimasaku, K., Malkan, M. A., et al. 2006, ApJ, 2, 7

[Kashino et al.(2019)Kashino, Lilly, Shibuya, Ouchi, &
  Kashikawa]Kashino2019arXiv
Kashino, D., Lilly, S. J., Shibuya, T., Ouchi, M., & Kashikawa, N. 2019,
  arXiv: 1909.09077

[Keating et al.(2019)Keating, Weinberger, Kulkarni, Haehnelt,
  Chardin, & Aubert]Keating2019
Keating, L. C., Weinberger, L. H., Kulkarni, G., et al. 2019, arXiv:
  1905.12640, 11, 1

[Kennicutt(1998)]Kennicutt1998
Kennicutt, R. C. 1998, STAR FORMATION IN GALAXIES ALONG THE HUBBLE SEQUENCE

[Konno et al.(2016)Konno, Ouchi, Nakajima, Duval, Kusakabe, Ono, &
  Shimasaku]Konno2016
Konno, A., Ouchi, M., Nakajima, K., et al. 2016, The Astrophysical Journal,
  823, 20

[Konno et al.(2018)Konno, Ouchi, Shibuya, Ono, Shimasaku, Taniguchi,
  Nagao, Kobayashi, Kajisawa, Kashikawa, Inoue, Oguri, Furusawa, Goto,
  Harikane, Higuchi, Komiyama, Kusakabe, Miyazaki, Nakajima, &
  Wang]Konno2018
Konno, A., Ouchi, M., Shibuya, T., et al. 2018, Publications of the
  Astronomical Society of Japan, 70

[Kulkarni et al.(2019)Kulkarni, Keating, Haehnelt, Bosman, Puchwein,
  Chardin, & Aubert]Kulkarni2019
Kulkarni, G., Keating, L. C., Haehnelt, M. G., et al. 2019, Monthly Notices
  of the Royal Astronomical Society: Letters, 485, L24

[Kusakabe et al.(2019)Kusakabe, Shimasaku, Momose, Ouchi, Nakajima,
  Hashimoto, Harikane, Silverman, & Capak]Kusakabe2019
Kusakabe, H., Shimasaku, K., Momose, R., et al. 2019, Publications of the
  Astronomical Society of Japan, 71, 1

[Laursen et al.(2019)Laursen, Sommer-Larsen, Milvang-Jensen, Fynbo,
  & Razoumov]Laursen2019
Laursen, P., Sommer-Larsen, J., Milvang-Jensen, B., Fynbo, J. P. U., &
  Razoumov, A. O. 2019, Astronomy & Astrophysics, 627, A84

[Leclercq et al.(2017)Leclercq, Bacon, Wisotzki, Mitchell, Garel,
  Verhamme, Blaizot, Hashimoto, Herenz, Conseil, Cantalupo, Inami, Contini,
  Richard, Maseda, Schaye, Marino, Akhlaghi, Brinchmann, &
  Carollo]Leclercq2017
Leclercq, F., Bacon, R., Wisotzki, L., et al. 2017, Astronomy &
  Astrophysics, 608, A8

[Livermore et al.(2017)Livermore, Finkelstein, &
  Lotz]Livermore2017
Livermore, R. C., Finkelstein, S. L., & Lotz, J. M. 2017, The Astrophysical
  Journal, 835, 113

[Madau(1995)]Madau1995
Madau, P. 1995, The Astrophysical Journal, 441, 18

[Madau & Haardt(2015)]Madau2015
Madau, P. & Haardt, F. 2015, Astrophysical Journal Letters, 813, L8

[Madau et al.(1998)Madau, Pozzetti, & Dickinson]Madau1998
Madau, P., Pozzetti, L., & Dickinson, M. 1998, The Astrophysical Journal, 498,
  106

[Malhotra & Rhoads(2004)]Malhotra2004
Malhotra, S. & Rhoads, J. E. 2004, The Astrophysical Journal, 617, L5

[Mason et al.(2019)Mason, Fontana, Treu, Schmidt, Hoag, Abramson,
  Amorin, Bradač, Guaita, Jones, Henry, Malkan, Pentericci, Trenti, &
  Vanzella]Mason2019
Mason, C. A., Fontana, A., Treu, T., et al. 2019, Monthly Notices of the
  Royal Astronomical Society, 485, 3947

[Mason et al.(2018)Mason, Treu, Dijkstra, Mesinger, Trenti,
  Pentericci, de Barros, & Vanzella]Mason2018
Mason, C. A., Treu, T., Dijkstra, M., et al. 2018, The Astrophysical Journal,
  856, 2

[Matsuoka et al.(2018)Matsuoka, Strauss, Kashikawa, Onoue, Iwasawa,
  Tang, Lee, Imanishi, Nagao, Akiyama, Asami, Bosch, Furusawa, Goto, Gunn,
  Harikane, Ikeda, Izumi, Kawaguchi, Kato, Kikuta, Kohno, Komiyama, Lupton,
  Minezaki, Miyazaki, Murayama, Niida, Nishizawa, Noboriguchi, Oguri, Ono,
  Ouchi, Price, Sameshima, Schulze, Shirakata, Silverman, Sugiyama, Tait,
  Takada, Takata, Tanaka, Toba, Utsumi, Wang, & Yamashita]Matsuoka2018
Matsuoka, Y., Strauss, M. A., Kashikawa, N., et al. 2018, The Astrophysical
  Journal, 869, 150

[McGreer et al.(2015)McGreer, Mesinger, & D'Odorico]McGreer2015
McGreer, I. D., Mesinger, A., & D'Odorico, V. 2015, Monthly Notices of the
  Royal Astronomical Society, 447, 499

[Momose et al.(2016)Momose, Ouchi, Nakajima, Ono, Shibuya,
  Shimasaku, Yuma, Mori, & Umemura]Momose2016
Momose, R., Ouchi, M., Nakajima, K., et al. 2016, Monthly Notices of the
  Royal Astronomical Society, 457, 2318

[Oke & Gunn(1983)]Oke1983
Oke, J. B. & Gunn, J. E. 1983, The Astrophysical Journal, 266, 713

[Ono et al.(2018)Ono, Ouchi, Harikane, Toshikawa, Rauch, Yuma,
  Sawicki, Shibuya, Shimasaku, Oguri, Willott, Akhlaghi, Akiyama, Coupon,
  Kashikawa, Komiyama, Konno, Lin, Matsuoka, Miyazaki, Nagao, Nakajima,
  Silverman, Tanaka, Taniguchi, & Wang]Ono2018
Ono, Y., Ouchi, M., Harikane, Y., et al. 2018, Publications of the
  Astronomical Society of Japan, 70, 2

[Ono et al.(2012)Ono, Ouchi, Mobasher, Dickinson, Penner, Shimasaku,
  Weiner, Kartaltepe, Nakajima, Nayyeri, Stern, Kashikawa, &
  Spinrad]Ono2012
Ono, Y., Ouchi, M., Mobasher, B., et al. 2012, The Astrophysical Journal,
  744, 83

[Ota et al.(2017)Ota, Iye, Kashikawa, Konno, Nakata, Totani,
  Kobayashi, Fudamoto, Seko, Toshikawa, Ichikawa, Shibuya, & Onoue]Ota2017
Ota, K., Iye, M., Kashikawa, N., et al. 2017, The Astrophysical Journal, 844,
  85

[Ota et al.(2012)Ota, Richard, Iye, Shibuya, Egami, &
  Kashikawa]Ota2012
Ota, K., Richard, J., Iye, M., et al. 2012, Monthly Notices of the Royal
  Astronomical Society, 423, 2829

[Ouchi et al.(2008)Ouchi, Shimasaku, Akiyama, Simpson, Saito, Ueda,
  Furusawa, Sekiguchi, Yamada, Kodama, Kashikawa, Okamura, Iye, Takata,
  Yoshida, & Yoshida]Ouchi2008
Ouchi, M., Shimasaku, K., Akiyama, M., et al. 2008, The Astrophysical Journal
  Supplement Series, 176, 301

[Ouchi et al.(2010)Ouchi, Shimasaku, Furusawa, Saito, Yoshida,
  Akiyama, Ono, Yamada, Ota, Kashikawa, Iye, Kodama, Okamura, Simpson, &
  Yoshida]Ouchi2010
Ouchi, M., Shimasaku, K., Furusawa, H., et al. 2010, The Astrophysical
  Journal, 723, 869

[Oyarzún et al.(2017)Oyarzún, Blanc, González,
  Mateo, & Bailey]Oyarzun2017
Oyarzún, G. A., Blanc, G. A., González, V., Mateo, M., & Bailey,
  J. I. 2017, The Astrophysical Journal, 843, 133

[Parsa et al.(2018)Parsa, Dunlop, & McLure]Parsa2018
Parsa, S., Dunlop, J. S., & McLure, R. J. 2018, Monthly Notices of the Royal
  Astronomical Society, 474, 2904

[Pentericci et al.(2011)Pentericci, Fontana, Vanzella, Castellano,
  Grazian, Dijkstra, Boutsia, Cristiani, Dickinson, Giallongo, Giavalisco,
  Maiolino, Moorwood, Paris, & Santini]Pentericci2011
Pentericci, L., Fontana, A., Vanzella, E., et al. 2011, The Astrophysical
  Journal, 743, 132

[Pentericci et al.(2018)Pentericci, Vanzella, Castellano, Fontana,
  De Barros, Grazian, Marchi, Bradac, Conselice, Cristiani, Dickinson,
  Finkelstein, Giallongo, Guaita, Koekemoer, Maiolino, Santini, &
  Tilvi]Pentericci2018
Pentericci, L., Vanzella, E., Castellano, M., et al. 2018, Astronomy &
  Astrophysics, 619, A147.

[Pentericci et al.(2014)Pentericci, Vanzella, Fontana, Castellano,
  Treu, Mesinger, Dijkstra, Grazian, Bradač, Conselice, Cristiani,
  Dunlop, Galametz, Giavalisco, Giallongo, Koekemoer, McLure, Maiolino, Paris,
  & Santini]Pentericci2014
Pentericci, L., Vanzella, E., Fontana, A., et al. 2014, The Astrophysical
  Journal, 793, 113

[Piqueras et al.(2017)Piqueras, Conseil, Shepherd, Bacon, Leclercq,
  & Richard]Piqueras2017
Piqueras, L., Conseil, S., Shepherd, M., et al. 2017, ArXiv, 1710.03554, 1

[Planck Collaboration et al.(2018)Planck Collaboration, Aghanim,
  Akrami, Ashdown, Aumont, Baccigalupi, Ballardini, Banday, Barreiro, Bartolo,
  Basak, Battye, Benabed, Bernard, Bersanelli, Bielewicz, Bock, Bond, Borrill,
  Bouchet, Boulanger, Bucher, Burigana, Butler, Calabrese, Cardoso, Carron,
  Challinor, Chiang, Chluba, Colombo, Combet, Contreras, Crill, Cuttaia,
  de Bernardis, de Zotti, Delabrouille, Delouis, Di Valentino, Diego,
  Doré, Douspis, Ducout, Dupac, Dusini, Efstathiou, Elsner, Enßlin,
  Eriksen, Fantaye, Farhang, Fergusson, Fernandez-Cobos, Finelli, Forastieri,
  Frailis, Franceschi, Frolov, Galeotta, Galli, Ganga, Génova-Santos,
  Gerbino, Ghosh, González-Nuevo, Górski, Gratton, Gruppuso,
  Gudmundsson, Hamann, Handley, Herranz, Hivon, Huang, Jaffe, Jones, Karakci,
  Keihänen, Keskitalo, Kiiveri, Kim, Kisner, Knox, Krachmalnicoff, Kunz,
  Kurki-Suonio, Lagache, Lamarre, Lasenby, Lattanzi, Lawrence, Jeune, Lemos,
  Lesgourgues, Levrier, Lewis, Liguori, Lilje, Lilley, Lindholm,
  López-Caniego, Lubin, Ma, Macías-Pérez, Maggio, Maino,
  Mandolesi, Mangilli, Marcos-Caballero, Maris, Martin, Martinelli,
  Martínez-González, Matarrese, Mauri, McEwen, Meinhold,
  Melchiorri, Mennella, Migliaccio, Millea, Mitra, Miville-Deschênes,
  Molinari, Montier, Morgante, Moss, Natoli, Nørgaard-Nielsen, Pagano,
  Paoletti, Partridge, Patanchon, Peiris, Perrotta, Pettorino, Piacentini,
  Polastri, Polenta, Puget, Rachen, Reinecke, Remazeilles, Renzi, Rocha,
  Rosset, Roudier, Rubiño-Martín, Ruiz-Granados, Salvati, Sandri,
  Savelainen, Scott, Shellard, Sirignano, Sirri, Spencer, Sunyaev, Suur-Uski,
  Tauber, Tavagnacco, Tenti, Toffolatti, Tomasi, Trombetti, Valenziano,
  Valiviita, Van Tent, Vibert, Vielva, Villa, Vittorio, Wandelt, Wehus, White,
  White, Zacchei, & Zonca]PlanckCollaboration6_2018arXiv
Planck Collaboration, Aghanim, N., Akrami, Y., et al. 2018, arXiv:,
  1807.06209, 1

[Price-Whelan et al.(2018)Price-Whelan, Sipőcz, Günther,
  Lim, Crawford, Conseil, Shupe, Craig, Dencheva, Ginsburg, VanderPlas,
  Bradley, Pérez-Suárez, de Val-Borro, Aldcroft, Cruz, Robitaille,
  Tollerud, Ardelean, Babej, Bach, Bachetti, Bakanov, Bamford, Barentsen,
  Barmby, Baumbach, Berry, Biscani, Boquien, Bostroem, Bouma, Brammer, Bray,
  Breytenbach, Buddelmeijer, Burke, Calderone, Rodríguez, Cara, Cardoso,
  Cheedella, Copin, Corrales, Crichton, D’Avella, Deil, Depagne, Dietrich,
  Donath, Droettboom, Earl, Erben, Fabbro, Ferreira, Finethy, Fox, Garrison,
  Gibbons, Goldstein, Gommers, Greco, Greenfield, Groener, Grollier, Hagen,
  Hirst, Homeier, Horton, Hosseinzadeh, Hu, Hunkeler, Ivezić, Jain,
  Jenness, Kanarek, Kendrew, Kern, Kerzendorf, Khvalko, King, Kirkby, Kulkarni,
  Kumar, Lee, Lenz, Littlefair, Ma, Macleod, Mastropietro, McCully, Montagnac,
  Morris, Mueller, Mumford, Muna, Murphy, Nelson, Nguyen, Ninan, Nöthe,
  Ogaz, Oh, Parejko, Parley, Pascual, Patil, Patil, Plunkett, Prochaska,
  Rastogi, Janga, Sabater, Sakurikar, Seifert, Sherbert, Sherwood-Taylor, Shih,
  Sick, Silbiger, Singanamalla, Singer, Sladen, Sooley, Sornarajah, Streicher,
  Teuben, Thomas, Tremblay, Turner, Terrón, Kerkwijk, de la Vega,
  Watkins, Weaver, Whitmore, Woillez, & Zabalza]TheAstropyCollaboration2018
Price-Whelan, A. M., Sipőcz, B. M., Günther, H. M., et al. 2018,
  The Astronomical Journal, 156, 123

[Rafelski et al.(2015)Rafelski, Teplitz, Gardner, Coe, Bond,
  Koekemoer, Grogin, Kurczynski, McGrath, Bourque, Atek, Brown, Colbert,
  Codoreanu, Ferguson, Finkelstein, Gawiser, Giavalisco, Gronwall, Hanish, Lee,
  Mehta, de Mello, Ravindranath, Ryan, Scarlata, Siana, Soto, &
  Voyer]Rafelski2015
Rafelski, M., Teplitz, H. I., Gardner, J. P., et al. 2015, The Astronomical
  Journal, 150, 31

[Richard et al.(2019)Richard, Bacon, Blaizot, Boissier, Boselli,
  NicolasBouché, Brinchmann, Castro, Ciesla, Crowther, Daddi, Dreizler,
  Duc, Elbaz, Epinat, Evans, Fossati, Fumagalli, Garcia, Garel, Hayes, Herrero,
  Humphrey, Jablonka, Kamann, Kaper, Kelz, Kneib, de Koter, Kudritzki, Langer,
  Lardo, Leclercq, Lennon, Mahler, Martins, Massey, Mitchell, Monreal-Ibero,
  Najarro, Opitom, Papaderos, Péroux, Revaz, Roth, Rousselot, Sander,
  Wagemann, Smail, Swinbank, Tramper, Urrutia, Verhamme, Vink, Walsh,
  Weilbacher, Wendt, Wisotzki, & Yang]Richard2019BlueMUSE
Richard, J., Bacon, R., Blaizot, J., et al. 2019, ArXiv:, 1906.01657

[Robertson et al.(2015)Robertson, Ellis, Furlanetto, &
  Dunlop]Robertson2015
Robertson, B. E., Ellis, R. S., Furlanetto, S. R., & Dunlop, J. S. 2015,
  Astrophysical Journal Letters, 802, L19

[Rykoff et al.(2015)Rykoff, Rozo, & Keisler]Rykoff2015arXiv
Rykoff, E. S., Rozo, E., & Keisler, R. 2015, arXiv.org, astro-ph.I, 1509.00870

[Salmon et al.(2015)Salmon, Papovich, Finkelstein, Tilvi, Finlator,
  Behroozi, Dahlen, Davé, Dekel, Dickinson, Ferguson, Giavalisco, Long,
  Lu, Mobasher, Reddy, Somerville, & Wechsler]Salmon2015
Salmon, B., Papovich, C., Finkelstein, S. L., et al. 2015, The Astrophysical
  Journal, 799, 183

[Santini et al.(2017)Santini, Fontana, Castellano, Criscienzo,
  Merlin, Amorin, Cullen, Daddi, Dickinson, Dunlop, Grazian, Lamastra, McLure,
  Michałowski, Pentericci, & Shu]Santini2017
Santini, P., Fontana, A., Castellano, M., et al. 2017, The Astrophysical
  Journal, 847, 76

[Santos et al.(2016)Santos, Sobral, & Matthee]Santos2016
Santos, S., Sobral, D., & Matthee, J. 2016, Monthly Notices of the Royal
  Astronomical Society, 463, 1678

[Schaerer et al.(2011a)Schaerer, De Barros, &
  Stark]Schaerer2011a
Schaerer, D., De Barros, S., & Stark, D. P. 2011a, Astronomy and
  Astrophysics, 536, 1

[Schaerer et al.(2011b)Schaerer, Hayes, Verhamme, &
  Teyssier]Schaerer2011
Schaerer, D., Hayes, M., Verhamme, A., & Teyssier, R. 2011b,
  Astronomy & Astrophysics, 531, A12

[Schenker et al.(2014)Schenker, Ellis, Konidaris, &
  Stark]Schenker2014
Schenker, M. A., Ellis, R. S., Konidaris, N. P., & Stark, D. P. 2014, The
  Astrophysical Journal, 795, 20

[Schmidt(1959)]Schmidt1959
Schmidt, M. 1959, The Astrophysical Journal, 129, 243

[Shapley et al.(2003)Shapley, Steidel, Pettini, &
  Adelberger]Shapley2003
Shapley, A. E., Steidel, C. C., Pettini, M., & Adelberger, K. L. 2003, The
  Astrophysical Journal, 588, 65

[Speagle et al.(2014)Speagle, Steinhardt, Capak, &
  Silverman]Speagle2014a
Speagle, J. S., Steinhardt, C. L., Capak, P. L., & Silverman, J. D. 2014, The
  Astrophysical Journal Supplement Series, 214, 15

[Stanway et al.(2008)Stanway, Bremer, & Lehnert]Stanway2008
Stanway, E. R., Bremer, M. N., & Lehnert, M. D. 2008, Monthly Notices of the
  Royal Astronomical Society, 385, 493

[Stark et al.(2009)Stark, Ellis, Bunker, Bundy, Targett, Benson, &
  Lacy]Stark2009
Stark, D. P., Ellis, R. S., Bunker, A., et al. 2009, The Astrophysical
  Journal, 697, 1493

[Stark et al.(2017)Stark, Ellis, Charlot, Chevallard, Tang, Belli,
  Zitrin, Mainali, Gutkin, Vidal-García, Bouwens, & Oesch]Stark2017
Stark, D. P., Ellis, R. S., Charlot, S., et al. 2017, Monthly Notices of the
  Royal Astronomical Society, 464, 469

[Stark et al.(2010)Stark, Ellis, Chiu, Ouchi, & Bunker]Stark2010
Stark, D. P., Ellis, R. S., Chiu, K., Ouchi, M., & Bunker, A. 2010, Monthly
  Notices of the Royal Astronomical Society, 408, 1628

[Stark et al.(2011)Stark, Ellis, & Ouchi]Stark2011
Stark, D. P., Ellis, R. S., & Ouchi, M. 2011, Astrophysical Journal Letters,
  728

[Steidel et al.(2000)Steidel, Adelberger, Shapley, Pettini,
  Dickinson, & Giavalisco]Steidel2000
Steidel, C. C., Adelberger, K. L., Shapley, A. E., et al. 2000, Astrophysical
  Journal Letters v.489, 532, 170

[The Astropy Collaboration et al.(2013)The Astropy
  Collaboration, Robitaille, Tollerud, Greenfield, Droettboom, Bray, Aldcroft,
  Davis, Ginsburg, Price-Whelan, Kerzendorf, Conley, Crighton, Barbary, Muna,
  Ferguson, Grollier, Parikh, Nair, Günther, Deil, Woillez, Conseil,
  Kramer, Turner, Singer, Fox, Weaver, Zabalza, Edwards, Azalee Bostroem,
  Burke, Casey, Crawford, Dencheva, Ely, Jenness, Labrie, Lim, Pierfederici,
  Pontzen, Ptak, Refsdal, Servillat, &
  Streicher]TheAstropyCollaboration2013
The Astropy Collaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013,
  Astronomy & Astrophysics, 558, A33

[Tilvi et al.(2014)Tilvi, Papovich, Finkelstein, Long, Song,
  Dickinson, Ferguson, Koekemoer, Giavalisco, & Mobasher]Tilvi2014
Tilvi, V., Papovich, C., Finkelstein, S. L., et al. 2014, Astrophysical
  Journal, 794, 1

[Totani et al.(2014)Totani, Aoki, Hattori, Kosugi, Niino, Hashimoto,
  Kawai, Ohta, Sakamoto, & Yamada]Totani2014
Totani, T., Aoki, K., Hattori, T., et al. 2014, Publications of the
  Astronomical Society of Japan, 66, 63

[Totani et al.(2006)Totani, Kawai, Kosugi, Aoki, Yamada, Iye, Ohta,
  & Hattori]Totani2006
Totani, T., Kawai, N., Kosugi, G., et al. 2006, Proceedings of the
  International Astronomical Union, 2, 265

[Treu et al.(2013)Treu, Schmidt, Trenti, Bradley, &
  Stiavelli]Treu2013
Treu, T., Schmidt, K. B., Trenti, M., Bradley, L. D., & Stiavelli, M. 2013,
  Astrophysical Journal Letters, 775, 5

[Ventou et al.(2017)Ventou, Contini, Bouché, Epinat,
  Brinchmann, Bacon, Inami, Lam, Drake, Garel, Michel-Dansac, Pello, Steinmetz,
  Weilbacher, Wisotzki, & Carollo]Ventou2017
Ventou, E., Contini, T., Bouché, N., et al. 2017, Astronomy &
  Astrophysics, 608, A9

[Verhamme et al.(2012)Verhamme, Dubois, Blaizot, Garel, Bacon,
  Devriendt, Guiderdoni, & Slyz]Verhamme2012
Verhamme, A., Dubois, Y., Blaizot, J., et al. 2012, A&A, 546, 111

[Verhamme et al.(2008)Verhamme, Schaerer, Atek, &
  Tapken]Verhamme2008
Verhamme, A., Schaerer, D., Atek, H., & Tapken, C. 2008, Astronomy &
  Astrophysics, 491, 89

[Verhamme et al.(2006)Verhamme, Schaerer, & Maselli]Verhamme2006
Verhamme, A., Schaerer, D., & Maselli, A. 2006, Astronomy and Astrophysics,
  460, 397

[Wilson(1927)]Wilson1927
Wilson, E. B. 1927, Journal of the American Statistical Association, 22, 209

[Wisotzki et al.(2016)Wisotzki, Bacon, Blaizot, Brinchmann, Herenz,
  Schaye, Bouché, Cantalupo, Contini, Carollo, Caruana, Courbot,
  Emsellem, Kamann, Kerutt, Leclercq, Lilly, Patrício, Sandin, Steinmetz,
  Straka, Urrutia, Verhamme, Weilbacher, & Wendt]Wisotzki2016
Wisotzki, L., Bacon, R., Blaizot, J., et al. 2016, A&A, 587, 98

[Yajima et al.(2012)Yajima, Li, Zhu, Abel, Gronwall, &
  Ciardullo]Yajima2012b
Yajima, H., Li, Y., Zhu, Q., et al. 2012, Astrophysical Journal, 754, 118

[Zheng et al.(2017)Zheng, Wang, Rhoads, Infante, Malhotra, Hu,
  Walker, Jiang, Jiang, Hibon, Gonzalez, Kong, Zheng, Galaz, &
  Barrientos]Zheng2017
Zheng, Z.-Y., Wang, J., Rhoads, J., et al. 2017, The Astrophysical Journal,
  842, L22



