© ESO 2021

1. Introduction

Over the last decade there have been increased efforts to characterize and study rest-frame ultraviolet (UV) emission lines from star-forming galaxies at increasingly higher redshifts, with a recent record holder at z = 11.09 presented by Jiang et al. (2020), as such studies provide extensive knowledge about the emitting galaxy population. Also at lower redshift and z ≈ 0, UV lines have recently provided the means necessary to leverage the extensive information provided by rest-frame optical emission, which is often unavailable at higher redshift, and they have provided analogs of high-redshift systems, being a key diagnostic for analyzing and understanding galaxy evolution in broader terms.

Rest-frame UV emission lines probe the physical conditions of the ionized gas, provide constraints on the physical properties of the emitting galaxies (such as from photo-ionization models), and have proven to be valuable probes of the surrounding environment of their host galaxies. Hence, they provide insights into general conditions for the galaxy formation and assembly at the targeted redshifts.

For instance, the [Si III] λ1883 + Si III] λ1892 and [C III] λ1907 + C III] λ1909 doublet components provide information about the electron density of the gas from which they are emitted (Osterbrock & Ferland 2006; Kewley et al. 2019a). The resonant C IVλ1548 + C IVλ1551 doublet is produced by highly ionizing radiation mainly from young stars, shocks, and/or active galactic nuclei (AGN), and it can be used to trace the physical conditions in the ionized gas and the interstellar medium (ISM) in star-forming galaxies as well as winds from the emitting massive O and B stars (Vidal-García et al. 2017; Byler et al. 2018; Berg et al. 2018; Feltre et al. 2020). By combining the [O III] λ1661 + O III] λ1666 and [C III] λ1907 + C III] λ1909 doublet fluxes, the C/O ratio in the galaxy can be approximated (Berg et al. 2018, 2016), and if information on the optical [O III] λ5007 flux can be obtained, the [O III] λ1661 + O III] λ1666 flux allows for a determination of the electron temperature (Berg et al. 2018), avoiding the need for detections of the faint auroral [O III] λ4363 emission (e.g., Sanders et al. 2016a, 2020; Bian et al. 2018). Determining the emission flux ratios of different UV line species and comparing them to grids of photoionization models predicting UV line fluxes provides estimates of, for instance, the metallicity and gas-phase abundances, the ionization parameter, and hydrogen number density (e.g., Gutkin et al. 2016; Feltre et al. 2016; Jaskot & Ravindranath 2016; Hirschmann et al. 2017, 2019; Byler et al. 2017, 2018, 2020; Nakajima et al. 2018; Plat et al. 2019; Kewley et al. 2019a). The He IIλ1640 emission likely includes both nebular emission and emission from stellar winds (Byler et al. 2018, 2020) and is therefore a valuable probe and diagnostic of these parameters, even though it has been proven challenging to reproduce the total observed He II emission from star-forming galaxies without invoking increased ionizing photon production from binary stars and/or X-ray binaries (e.g., Steidel et al. 2016, 2018; Berg et al. 2019a; Nanayakkara et al. 2019; Senchyna et al. 2021).

The intrinsically strongest UV emission line H Iλ1216 (Lyα) is the most common spectroscopic probe of high-z galaxies. But, due to the high Lyα absorption cross-section in neutral hydrogen, Lyα photons scatter resonantly in the ISM and circumgalactic medium (CGM). This scattering results in Lyα emission being more extended than the stellar continuum emission (Steidel et al. 2011; Momose et al. 2014; Wisotzki et al. 2016, 2018; Leclercq et al. 2017, 2020) and imprints kinematic and spatial properties of the ISM and CGM gas distribution into the observable Lyα signal from the Lyα emitters (LAEs; used here for any galaxy with detected Lyα flux irrespective of the rest-frame equivalent width, EW₀). The scattering furthermore makes Lyα photons more prone to dust extinction, which can lead to a complete suppression of Lyα from galaxies (e.g., Dijkstra et al. 2011; Laursen et al. 2011, 2019; Dijkstra 2017; Kusakabe et al. 2020, see also Sect. 7). By combining the Lyα observables with inferences from the other rest-frame UV lines we may constrain this interplay between the galaxies gas and the Lyα radiation field. This is of particular value for predicting Lyα observables from other UV lines during the Epoch of Reionization (EoR) at z ≳ 6, where the neutral hydrogen in the intergalactic medium (IGM) scatters Lyα out of the observers line of sight (e.g., Treu et al. 2013; Pentericci et al. 2014; Tilvi et al. 2014; de Barros et al. 2017; Caruana et al. 2018; Kusakabe et al. 2020). The differences between observed and expected Lyα can therefore be used to constrain the temporal evolution of the degree of ionization in the universe during the EoR.

Hence, in addition to providing detailed knowledge about the detected galaxies themselves, rest-frame UV emission lines other than Lyα can serve as redshift confirmation of galaxies in neutral environments in the EoR when Lyα is completely absent and provide probes of the systemic redshift of LAEs when Lyα is offset due to to the scattering in the ISM and CGM. To assess and quantify the feasibility of such approaches, the characteristics of these secondary UV emission lines need to be explored and understood for LAEs at non-EoR redshifts, that is z < 6. Such efforts to characterize and study UV emission lines at both intermediate and high redshift will support the analysis of data from upcoming space-based sensitive infrared observatories like the James Webb Space Telescope (JWST) and the Nancy Grace Roman Space Telescope (formerly known as WFIRST) which are expected to revolutionize studies of galaxy evolution and assembly deep into the epoch of reionization (EoR).

Previously, sources with UV emission lines have been collected through three main methods (see also Byler et al. 2020). First, creating composite spectra from stacking large samples of individual galaxy spectra provided some of the first results (e.g., Shapley et al. 2003) and has since been expanded (Steidel et al. 2016; Feltre et al. 2020). Secondly, objects with extreme emission line ratios and fluxes high enough for direct detection have also been found (e.g., Erb et al. 2010; Stark et al. 2014; Amorín et al. 2017). Thirdly, targeting gravitationally lensed sources uses nature’s own telescopes to improve the observation’s sensitivity to be able to detect the often intrinsically faint UV emission at high redshift (e.g., Bayliss et al. 2014; Patrício et al. 2016; Schmidt et al. 2017; Rigby et al. 2018a,b; Berg et al. 2018). To gain insight in the physical properties governing galaxy formation from the EoR to the main epoch of star formation at z ≈ 2, in this paper we used a fourth method. We obtain samples of rest-frame UV emitters, by exploring observations from the Multi Unit Spectroscopic Explorer (MUSE; Bacon et al. 2010) at ESO’s Very Large Telescope (VLT). A similar approach, that is searching for UV emission in large parent samples with rest-frame UV spectroscopy has also been pursued by, for example, Maseda et al. (2017) and Nanayakkara et al. (2019) using MUSE-Deep data, Saxena et al. (2020) using the VANDELS data, and Nakajima et al. (2018) and Le Fevre et al. (2019) using the VIMOS Ultra-Deep Survey. Basing these studies on samples of spectroscopically confirmed objects ensures confident redshifts aiding the efficiency of the search for the intrinsically faint UV emission. Additionally, studying objects from samples based on emission line objects as opposed to objects with redshift determined from continuum spectroscopy, allows probing populations of fainter galaxies and ensures the presence of recent star formation responsible for the nebular UV line emission.

Here we explore data taken with MUSE of 2052 spectroscopically confirmed emission line galaxies at 1.5 < z < 6.4 including all LAEs and potential [C III] λ1907 + C III] λ1909 emitters. The MUSE data are presented in Sect. 2. Sections 3 and 4 explain how we selected and determined the redshifts of our parent sample and how we extracted optimal one-dimensional (1D) spectra of each of these sources. Section 5 describes our emission line template fitting approach (detailed in Appendix A) used to search for UV emission red-wards of Lyα in the 1D spectra. This results in a large sample of UV line emitters of both LAEs and non-LAEs which we explore the physical parameters of in the remainder of the paper. We infer EW₀ of the emission lines (Sect. 6) and investigate correlations between them and the emission line fluxes with properties of the LAEs in our sample (Sect. 7). We explore the Lyα emission line velocity offsets with respect to the secondary UV emission lines which provide systemic redshifts (Sect. 8), and assess the range of electron densities for objects showing [Si III] λ1883 + Si III] λ1892 and [C III] λ1907 + C III] λ1909 (Sect. 9). In Sect. 10 we determine the gas-phase abundances of systems with multiple UV emission features detected. Finally, we explore how further insight can be gained on the physical parameters of the targeted galaxies from photoionization models by introducing “PhotoIonization Model Probability Density Functions” (PIM-PDFs) in Sect. 11, before summarizing and concluding our study in Sect. 12. With the paper, we provide a catalog of the complete set of UV emission line measurements from the 2052 galaxies studied (described in Appendix B). Throughout the paper we compare these measurements to a collection of measurements from the literature which we also provide a catalog of Appendix C.

In the remainder of this paper, we use the short notations for the rest-frame UV emission lines red-wards of Lyα listed in Table 1. The minimum required energy of ionizing photons responsible for these UV lines are also quoted together with the expected size of the emission doublet line flux ratios. We use AB magnitudes (Oke & Gunn 1983) and assume flat cosmological parameters of H₀ = 70 km s⁻¹ Mpc⁻¹, Ω_m = 0.3, and Ω_Λ = 0.7.

2. Data

In this paper we search for UV emission lines in a compilation of the MUSE consortium guaranteed time observations (GTO) data. The MUSE integral field spectrograph provides spatial and spectral coverage of the 1 arcmin² field-of-view (FoV). The spatial sampling in the reduced data cubes is . Each of these spaxels contains a spectrum spanning the optical wavelength range from 4800 Å to 9300 Å at a wavelength sampling of 1.25 Å per volume element, that is per voxel. The spectral resolution element of the MUSE data is roughly 2.5 Å across the full wavelength range such that the spectral resolution R = λ/Δλ ≈ 3000 (R ≈ 1800 in the blue up to 4000 in the red, Weilbacher et al. 2020). The data sets explored in this study come from the MUSE GTO team’s “wedding cake approach” of a shallow, medium-deep and deep survey of multiple MUSE pointings. We describe each of these data sets in the following subsections. The general outline and the location of each of the MUSE pointings are shown in Fig. 1.

Fig. 1. Location of the analyzed MUSE pointings in the GOODS-South (top), HUDF (center) and COSMOS (bottom) regions. The MUSE-Wide pointings are shown in orange, the UDF mosaic in blue and the UDF10 pointing in red. The location of the 2052 sources searched for rest-frame UV emission lines are marked with green circles. The gray scale images are the 3D-HST combined near-infrared detection images from https://3dhst.research.yale.edu/ (top and bottom) and the HLF F775W version 2.0 image from https://archive.stsci.edu/prepds/hlf/ (center).

2.1. MUSE-Wide

The wide and shallow component of the GTO wedding cake approach is the MUSE-Wide survey. It consists of 100 1 × 1 arcmin² MUSE pointings distributed over the COSMOS (Scoville et al. 2007) and GOODS-South (Dickinson et al. 2003; Giavalisco et al. 2004) CANDELS footprints (Grogin et al. 2011; Koekemoer et al. 2011).

The first catalog and data release of the MUSE-Wide survey were presented by Herenz et al. (2017) and Urrutia et al. (2019), respectively. These studies analyze 44 MUSE pointings of the MUSE-Wide data over the GOODS-South footprint. Here we include data from these 44 fields adding the remaining 47 MUSE pointings of the full MUSE-Wide data set. Formally, the MUSE-Wide survey consists of 100 MUSE pointings. But when referring to MUSE-Wide throughout this paper, we are excluding the nine pointings in the Hubble Ultra Deep Field (HUDF; Beckwith et al. 2006), as these fields are included in the deeper MUSE UDF mosaic described in Sect. 2.2. The MUSE-Wide data were obtained through a series of MUSE GTO observing runs from September 2014 through September 2017. Each of the 91 MUSE-Wide pointings have a depth of 1 h and were observed on nights with a seeing of roughly . Details of the seeing properties for the first data release from MUSE-Wide are given in Sect. 3.2.5 of Urrutia et al. (2019). The remaining MUSE-Wide fields were observed under similar conditions.

The MUSE-Wide data were reduced using version 1.0 (or an early development version) of the data processing pipeline (DPP) for the MUSE instrument (Weilbacher et al. 2014, 2020) released in December 2014 to support the first observing runs with the seeing-limited wide-field mode of MUSE. The DPP was complemented by enhanced so-called slice-based sky-subtraction, a self-calibration of the initial reductions using the self-calibration procedure from the MUSE Python Data Analysis Framework (MPDAF; Conseil et al. 2016; Piqueras et al. 2017), and accounting for the effect of varying shape of the LSF on the sky residuals with the Zurich Atmosphere Purge (ZAP version 1.0, Soto et al. 2016). Even though newer versions of the DPP have been released since 2014, all 91 MUSE-Wide pointings were reduced using the same versions of the DPP and the complementing software for consistency. For details on the MUSE-Wide data reduction see Urrutia et al. (2019). The reduced MUSE-Wide data on average reach a 1σ emission line detection sensitivity of 1 × 10⁻¹⁹ erg s⁻¹ cm⁻² Å⁻¹ as estimated from point source and LAE insertion experiments (Herenz et al. 2019) assuming a 1 arcsec² aperture.

2.2. MUSE-Deep UDF mosaic

Complementing the wide and relatively shallow data from MUSE-Wide, the GTO wedding cake program contains a medium depth survey of nine pointings at ten hours depth over the HUDF which we refer to as the UDF mosaic (blue squares in the top and central panels of Fig. 1). The data over the UDF analyzed in the current study were first presented by Bacon et al. (2017) and Inami et al. (2017). As described in Sect. 3 below, we perform an independent and self-consistent source selection in all data analyzed in our study. Hence, we are relying on these source lists instead of the source catalog presented by Inami et al. (2017). The UDF mosaic data were taken through a series of MUSE GTO observing runs from September 2014 through February 2016. Our source selection is based on an updated and improved reduction of these UDF data with the MUSE DPP which will be presented by Bacon et al. (in prep.). Like the MUSE-Wide data reduction, the basic DPP reduction was complemented by the self-calibration and improved sky subtraction from MPDAF and ZAP for the reduction of the deeper data. For details on the UDF mosaic (and UDF10) data reduction see Bacon et al. (in prep.). As shown by Bacon et al. (2017) the first version of the UDF mosaic data reach a 1σ emission line flux sensitivity of 5.5 × 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹ for a 1 arcsec² aperture averaged over the most sensitive part of the MUSE throughput curve at 7000−8500 Å.

2.3. MUSE-Deep UDF10

The final layer of the MUSE GTO wedding cake program is the field referred to as UDF10. UDF10 is a single 1 × 1 arcmin² MUSE pointing at 31 h depth within the HUDF (red square in top and central panels of Fig. 1). The UDF10 data are also part of the updated release of the MUSE-Deep data presented by Bacon et al. (in prep.) and were reduced using the same setup as the UDF mosaic data. For the first version of the MUSE-Deep data presented by Bacon et al. (2017), the estimated 1σ emission line flux for the UDF10 data is 2.8 × 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹.

3. Object selection in the MUSE data

To establish a sample of line emitters in the reduced MUSE GTO data cubes, all 101 individual data cubes were searched for emission lines using the dedicated Line Source Detection and Cataloguing (LSDCat) software described by Herenz & Wisotzki (2017). Apart from the flux cube LSDCat requires a variance cube for the data too. As the variances obtained by formal error propagation in the MUSE DPP underestimate the true uncertainties of the data due to the voxel resampling in the cube construction, we estimated empirically calibrated effective variances for all data cubes. In short, the effective variance cubes are constructed by first (i) measuring the typical variances between individual sky-voxels to capture the small-scale systematics like imperfect sky subtraction and flat fielding, then (ii) account for the covariances from the noise propagation in the DPP by rescaling the noise level by the average covariance level estimated by pulling a random noise cube (Gaussian with mean 0 and variance 1) through the DPP, and lastly (iii) assuming that the effective noise is constant (modulo the number of exposures) across each individual wavelength slice of the data cubes, that is assuming background-limited noise. For further details on the effective variance procedure see Sect. 3.2.4 of Urrutia et al. (2019) and Weilbacher et al. (2020). Furthermore, LSDCat requires continuum subtracted data cubes for line detection. These cubes were generated from the reduced data cubes by subtracting a 151 pixels (≈189 Å) wide running median from each of the spaxels.

The untargeted search for emission lines with LSDCat is based on these continuum subtracted versions of the reduced data cubes and the effective variance cubes. LSDCat then performs a 3D template match over all voxels in the data cube. To optimize the recovery of Lyα emitters, we used a search template with a Gaussian spectral component of fixed full width at half maximum (FWHM) set to 250 km s⁻¹ and a two-dimensional (2D) spatial Gaussian component with FWHM equal to the width of the point spread function (PSF) of the observations in the individual data cubes. Hence, the detection significance is slightly biased towards narrow emission lines. However, this does not affect our independent search for rest-frame UV lines (Sect. 5) in the sources recovered by LSDCat. The untargeted emission line detections from LSDCat with a signal-to-noise ratio (S/N) above five (6.4 for the first 24 MUSE-Wide fields, see Urrutia et al. 2019) were grouped and then classified with the visual inspection tool QtClassify (Kerutt 2017; Herenz et al. 2017). Herenz & Wisotzki (2017) provide the analytic approximation for the minimum line flux recoverable by LSDCat at a given detection threshold¹

(1)

Considering S/N = 5 detections and inserting the Δλ = 1.25 Å wavelength sampling, σ_z = 1.46 pixels at 5000 Å corresponding to the spectral detection template FWHM of 250 km s⁻¹, and a spatial source model width of σ_G = 1.84 pixels corresponding to a spatial Gaussian width of the flux sensitivities quoted in Sect. 2 translate into 5σ limiting line fluxes of roughly 1.3 × 10⁻¹⁷, 7.2 × 10⁻¹⁸, and 3.7 × 10⁻¹⁸ erg s⁻¹ cm⁻², for the MUSE-Wide, UDF mosaic, and UDF-10 data, respectively.

For the source classification, two (or sometimes three) people independently categorized all line-sets found by LSDCat by matching the shape and spectral location with rest-frame emission line lists. The independent classifications were then consolidated in plenum with a third independent consolidator ensuring at least three independent assessments of all detected emission lines from the untargeted search. The leading line of each line emitter, that is the associated line with the highest 3D LSDCat S/N was assigned a type (for instance, Lyα, [O II] or Hα) and a confidence from C = 0 to C = 3. A confidence of 1 refers to line emitters with a single trustworthy line detection of uncertain nature, C = 2 refers to a classification with high confidence, where there is support from line profiles (for example, a skewed asymmetric profile for Lyα emitters) or secondary S/N < 5 lines. For objects with C = 3 there is little room for doubt about the line classification, either due to characteristic line profiles or multiple line detections at S/N > 5. A confidence of 0 was given on rare occasions where it was hard to tell if the emission feature was real or not.

The QtClassify classification of the detected line emitters from LSDCat from all 101 MUSE data cubes leads to the self-consistent parent line emitter catalog from MUSE-Wide, the UDF mosaic and UDF10 used here. This catalog therefore differs from both the line emitters presented by Inami et al. (2017) and the source catalog based on the improved rereduction of the MUSE-Deep data by Bacon et al. (in prep.). The Inami et al. (2017) emission line sources were detected using the ORIGIN software (Mary et al. 2020) in the data cubes from the first MUSE-Deep data release (Bacon et al. 2017). The updated release of the UDF mosaic and UDF10 catalogs by Bacon et al. (in prep.) merges the independent LSDCat and ORIGIN emission line source detections performed on the improved reductions of the MUSE-Deep data complemented by even deeper MUSE data in the HUDF from the MUSE eXtremely Deep Field (MXDF; see Bacon et al. 2021). Bacon et al. (in prep.) also present continuum-selected non-emission line sources and sources at z < 1.5. The self-consistent LSDCat-based catalog of line emitters at z > 1.5 presented as part of our work is well suited for statistical studies. However, in terms of source completeness the MUSE-Deep catalog from the UDF mosaic and UDF10 by Bacon et al. (in prep.) will supersede it.

Each individual LSDCat detection in MUSE-Wide, the UDF mosaic or UDF10 was assigned a unique ID which we use throughout this paper when referring to and discussing individual objects. For the sources from the MUSE-Wide pointings we follow the ID structure of the first data release of MUSE-Wide (Urrutia et al. 2019) where each ID contains nine digits using the format “ABBCCCDDD”. Here A refers to the region in which the object was detected (1 = GOODS-South, 2 = COSMOS, 3 = northern HUDF parallel, and 4 = southern HUDF parallel), BB refers to the arbitrary numbering of the individual MUSE pointings, CCC refers to the LSDCat object identifier for that particular data cube, and DDD refers to a running number for the detected lead line of the source. The sources in the UDF mosaic have a nine digit ID using the format “ACCCCDDDD”. Here A = 6 refers to the UDF mosaic, CCCC is the four-digit LSDCat object identifier, and DDDD is the running number for the detected lead line of the source. Similarly, the UDF10 objects were assigned IDs using the format “A2CCCDDDD” where A = 7 refers to UDF10, and CCC and DDDD again indicate the LSDCat object identifier and the running number of the lead line.

To focus on potential rest-frame UV line emitters and for completeness, we selected all objects with a lead line at a redshift above 1.5 from the parent catalog irrespective of the confidence assignment. This includes all objects with Lyα as the lead line corresponding to z ≳ 2.9 for the MUSE wavelength coverage and all objects in the so-called MUSE redshift desert at z ≈ 1.5−2.9. At z ≈ 1.5 C III enters the MUSE wavelength range, whereas the prominent optical [O II] λ3726 + [O II] λ3729 doublet moves out of it. Hence, neither Lyα nor [O II] λ3726 + [O II] λ3729 is available for redshift determination in the redshift desert. In total there are 2197 unique IDs in the 101 MUSE pointings satisfying these selection criteria (1120 from MUSE-Wide, 842 from the UDF mosaic and 235 from UDF10). The location of these emission line sources are marked by the green dots in Fig. 1 and do not account for the overlaps between individual pointings.

4. Optimal spectral extraction

Our search for faint UV emission line features in the parent sample of the selected line emitters from the MUSE GTO data is performed on 1D spectra extracted from the 3D data cubes. Here we use optimally extracted 1D spectra obtained with the dedicated tool for Three-Dimensional Optimal Spectral Extraction (TDOSE²; Schmidt et al. 2019). In short, TDOSE bases the spectral extraction on a 2D spatial model of each object in the FoV. This model is based on higher-resolution imaging from, for instance, the Hubble Space Telescope (HST). Hence, it is implicitly assumed that any emission line fluxes extracted from these spectra are well represented by the extent of the continuum emission in the modeled HST images, which is a poor assumption for extended emission like Lyα but reasonable for UV emission lines. Using the morphological HST models as object templates, TDOSE performs a simultaneous 3D template match of all sources in the FoV, scaling source fluxes at each wavelength layer in the data cube accounting for contaminating flux by neighboring sources. Hence, TDOSE is based on an approach similar to PampelMuse (Kamann et al. 2013; Kamann 2018), except that TDOSE focuses on extended sources (galaxies), whereas PampelMuse was developed for spectral extraction of point-sources (stars).

All objects in our parent sample have spectra extracted using a single multivariate Gaussian source model for both the object of interest and each of the contaminants in the FoV. For source extractions from the 91 MUSE-Wide data cubes the contamination models were based on the Skelton et al. (2014) 3D-HST photometric catalog, as this catalog provides self-consistent source catalogs over GOODS-South, the HST UDF parallel fields and the COSMOS region mapped by the MUSE-Wide survey (see Fig. 1). In GOODS-South and COSMOS we modeled the source morphology in the available CANDELS HST F814W images corresponding to the Skelton et al. (2014) photometry. For the two HUDF parallel fields we based the models on the HST F160W images to ensure full coverage of the 2 × 2 MUSE pointings using the corresponding CANDELS and 3D-HST imaging and to ensure access to morphological priors from the Skelton et al. (2014) catalog for the modeling of the HST images. For extractions in the HUDF region we based the spatial models on the source catalog by Rafelski et al. (2015) and the version 1.5 HST F775W images from the Hubble Legacy Field (HLF³; Illingworth et al. 2016). No additional exposures were added to the F775W images for the recent HLF version 2.0 images presented by Whitaker et al. (2019). Spectra for sources with no detections in the ancillary imaging, that is where either faint or no photometric counterparts could be identified, were extracted using a morphological model identical to a point-source convolved with the PSF of the observations.

After an initial run of fully automatic spectral extractions, the resulting 1D spectra were visually inspected to identify suboptimal source models or cases where the source flux scalings were unreliable. Roughly 2% (42/2197) of the spectra were selected for reextraction with more careful attention to source location, source numbers, FoV extent, etc. Among the 2197 initial spectra only five were of LAEs with z > 6.44. At these redshifts all of the considered secondary emission lines are redshifted beyond the red cutoff of the MUSE instrument at λ ≈ 9300 Å. These five spectra were therefore not searched for UV emission lines. For four objects (126042110, 602121764, 609223654, and 613534254) we rely on aperture spectra, as no satisfactory spatial morphological multivariate Gaussian or PSF-based extraction could be obtained when including the complete contamination model. They were also obtained with TDOSE, using an extraction aperture with a radius of , which corresponds roughly to the FWHM of the MUSE PSF in the deeper MUSE data (Bacon et al. 2017). Two objects (158002004, 601931670) were located directly behind bright foreground objects hampering a satisfactory spectral extraction. These two objects were also excluded in the further analysis. Finally, two objects (208014258 and 600341002) were removed from the sample due to severe contamination and unreliable flux scalings caused by neighboring stars.

As can be seen in Fig. 1, the UDF mosaic (blue pointings) overlaps with parts of the MUSE-Wide coverage (orange pointings) and the UDF10 field (red square) is fully within the UDF mosaic. Hence, the UDF mosaic and UDF10 contain duplicates of sources in MUSE-Wide and the UDF mosaic, respectively. These duplicates were identified by searching for multiple objects within a search radius of of all positions in the main catalog with lines identified at similar wavelengths. The search radius corresponds to half the approximate coordinate precision between the MUSE detections and the HST reference images of . We excluded the spectra of the UDF mosaic (MUSE-Wide) sources when there was a UDF10 (UDF mosaic) counterpart within this search radius at the same redshift, which corresponds to a total of 120 (16) sources. Hence, the final sample of emission line sources with 1.5 < z < 6.4 from MUSE-Wide, the UDF mosaic and UDF10 analyzed in the remainder of this paper amounts to 2052 unique objects of which 1997 are LAEs. The thick green histogram in Fig. 2 shows the redshift distribution of these objects compared to the independent MUSE-Deep UDF mosaic and UDF10 catalog presented by Inami et al. (2017, thin blue histogram) and the MUSE-Wide DR1 catalog by Urrutia et al. (2019, thin orange histogram). The set of objects at z > 1.5 from the latter is a subsample of the parent sample studied here. The lack of Lyα and [O II] λ3726 + [O II] λ3729 for redshift identification resulting in the MUSE redshift desert is clearly visible and marked by the gray shaded region. Figure 3 presents the TDOSE spectra of two example objects. Further examples of objects representing the breadth of data and supporting discussions of individual objects in the remainder of the paper are shown in Appendix E.

Fig. 2. Redshift distributions of the sources studied here (thick green histogram), the subsample of them with UV emission lines detected by FELIS (filled dark green histogram), the comparison sample of UV line emitters from the literature described in Appendix C (dotted black histogram), the MUSE-Wide DR1 catalog by Urrutia et al. (2019, thin orange histogram), and the independent MUSE-Deep UDF mosaic and UDF10 catalog presented by Inami et al. (2017, thin blue histogram). The MUSE redshift desert where neither [O II] λ3726 + [O II] λ3729 (0.3 < z < 1.5) nor Lyα (z > 2.9) emission is available for redshift identification in the MUSE wavelength range is marked by the gray band.

Fig. 3. Examples of TDOSE extracted 1D MUSE spectra (blue) of the two LAEs 201004081 (z = 3.08; top) and 602922055 (z = 3.44; bottom) with their 1σ error spectrum indicated by the blue-shaded region. Further examples can be found in Appendix E. Each of the smaller panels shows a UV emission line region as marked by the black boxes in the bottom panels, which show the full MUSE wavelength range. Vertical gray lines mark the approximate location of emission features based on the redshift of the LSDCat lead line. No continuum subtraction has been performed on the spectra. The UV emission lines detected by the software FELIS (Sect. 5) are marked by the red curves. The solid red curves show FELIS detections with S/N > 3 and a velocity offset below 1000 km s−1 deemed trustworthy in the visual inspection process. Dashed red curves show less secure detections, which were however deemed reliable. Dotted red curves show S/N > 3 and |Δv|< 1000 km s−1 FELIS detections classified as unreliable. In each panel 2× the median data cube noise in the respective fields (MUSE-Wide, MUSE mosaic and UDF10) are shown for reference as the black filled curves offset to the bottom of each panel for clarity. These median noise spectra are also shown in Fig. A.2.

5. Searching for UV emission lines in 1D spectra

Knowing the (approximate) redshift for each of the 2052 unique objects, we searched the TDOSE spectra for signs of rest-frame UV emission lines through template matching using the publicly available Python software for Finding Emission Lines in Spectra (FELIS⁴ described in Appendix A). In short, FELIS matches a set of predefined emission line (doublet) templates independently to each spectrum via cross-correlation and χ² minimization. By cross-correlating each template (𝒯) with the spectrum around the expected location of the UV emission line, the template flux scaling that minimizes the disagreement with the data is estimated (α in Eq. (A.4)). Calculating the S/N of the best match of the cross-correlations provides an estimate for (S/N)_𝒯, max for each template and its flux scaling. The matched emission line (doublet) template with the highest (S/N)_max then provides the best match to the considered part of the observed spectrum overall. We refer to this best-match S/N as the “FELIS S/N” of the detected emission lines in the remainder of this work. Here we focus our search on the rest-frame UV emission lines red-wards of Lyα listed in Table 2.

All emission line template matches were performed independently allowing for individual velocity offsets with respect to the LSDCat lead line of the spectra, which provides the selection redshift as described in Sect. 3. The lead line redshifts (and the targeted UV lines) are not guaranteed to be at systemic redshift, especially for objects where the lead line is the resonant C IV or Lyα line, which can be offset from systemic by hundreds of km s⁻¹, as indicated by the location of the vertical gray lines based on the LSDCat lead line redshift shown in Fig. 3. Each spectrum was therefore searched around the expected location of the central UV line wavelengths (Table 2) in the range [4790/(1+z) + 10 Å, 9310/(1+z) − 10 Å]. Hence, if a central doublet wavelength is outside this range it will not be recovered by our search. This means that we search for the different rest-frame UV emission lines in the effective redshift ranges provided in Table 2.

The right part of Table 2 lists the characteristics of the emission line templates used to search for each of the emission features in the 1D TDOSE spectra. FELIS provides tools for generating search templates as described in Appendix A.1. Here we generated Gaussian templates for all lines with a given width (σ_Gauss in steps of Δ as given in Table 2). For the emission line doublets we generated templates with two Gaussian components fixing the separation between them. However, the doublet flux ratios were free to vary as listed in Table 2. By using a fixed set of Gaussian templates we obtain robustness and efficiency when searching for undiscovered UV emission in the more than two thousand spectra compared to performing parametric fits of Gaussians at all possible locations in the spectra.

Of the 2052 objects searched for UV emission lines we obtained an initial candidate list of 705 objects with at least one potential UV line detection with FELIS S/N > 3 and a velocity offset with respect to the catalog redshift of less than 1000 km s⁻¹. Despite including the effective noise error information in the FELIS matches, a large fraction of these potential detections were spurious emission lines, where the cross-correlation latched on to sky-line residuals or minor variations in the overall spectrum matching template doublet spacings. Using FELIS S/N > 5 reduces the number of objects with at least one candidate UV line detection to 167 objects removing the majority of the spurious detections but also discarding reliable detections. To account for the spurious detections while still recovering reliable detections at 3 < S/N < 5, we therefore visually vetted all 705 > S/N > 3 emitter candidates to determine the amount and type of reliable FELIS template fits. Figure 3 shows examples of UV emission line detections deemed reliable (red solid curves over-plotted on the blue TDOSE spectra) and FELIS detections discarded in the vetting process (dashed and dotted red curves). Hence, after visually vetting the 705 (167) objects with at least one potential S/N > 3 (5) UV emission line detection we deemed 103 (57) of these objects to have reliable FELIS template matches to the MUSE TDOSE spectra. A summary of the 103/705 objects with reliable UV emission line detections is presented in Table 3 and their redshift distribution is shown as the filled dark green histogram in Fig. 2. In Table 3 we list the number of detected lines, the number of available objects to search with FELIS and the corresponding fraction of detections for both the full sample (top), objects in the redshift desert (central), and the LAE sample (bottom) for all 101 fields, the MUSE-Wide pointings, the UDF mosaic and the deep UDF10 pointing. Of the listed 103 objects with UV line detections, 71 and 80 lines have S/N > 5 and 3 < S/N < 5, respectively. The objects with line detections in the redshift desert generally have high-confidence redshift determinations as 84% of the sources (41/49) have at least one detection with S/N > 5. Furthermore, 53% of the objects (26/49) have multiple lines detected, where at least one line has S/N > 5. Figure E.5 shows two example of such redshift desert objects with multiple UV lines detected at high confidence. Appendix B and Table B.1 describe the value-added catalog of the full sample of objects searched for UV emission lines including the corresponding fluxes, flux ratios, and EWs provided with this paper. This catalog also includes the relevant measurements of the Lyα properties determined by Kerutt et al. (2021) described in Sect. 7.

Figure 4 shows the UV emission line detection fractions for each of the subsamples. In general, there is an increase of the fraction of objects with detected UV emission lines red-wards of Lyα as the depth of the data is increased (1, 10, and 31 h depth for MUSE-Wide, UDF mosaic and UDF10, respectively). In some cases the low number of detections might be responsible for the potential trends as indicated by the error bars showing the 95% confidence intervals for the Clopper–Pearson interval (Clopper & Pearson 1934). These intervals are conservative as they are extracted based directly on the cumulative probability function of the binomial distribution but reflect the uncertainty of each fraction. An increase in the fraction of faint emission lines as a function of survey depth will only occur if there is a relative change in the shape of the object’s (emission line) luminosity function for the different samples. For example, if the truncated LAE luminosity function shape (slope) was the same as that of the UV emitter luminosity function at the various depths, the fraction of objects with detected UV emission lines above the 3σ threshold should remain the same irrespective of survey depth. Given the limited sample size and survey depth Stroe et al. (2017) fixed the faint-end slope for their C III LF fit suggesting that it is close to a scaled-down version of the Lyα LF. This is not what we see in the MUSE data, implying that the shape of the luminosity functions of the subsamples of objects emitting the detected UV emission lines differ in shape compared to the luminosity function of the parent (LAE) sample. A steeper faint-end slope of the UV emission line luminosity function could replicate the observed increase in the fraction of UV emitters with increasing depth.

Fig. 4. Percentage of objects with detected secondary UV emission lines red-wards of Lyα for the full sample (circles) and the subsample of LAEs (squares). The percentages are provided for all fields studied, the MUSE-Wide data (1 h depth), the UDF mosaic data (10 h depth) and the UDF10 pointing (31 h depth). The error bars indicate the 95% confidence Clopper–Pearson interval of the binomial distribution for k UV line detections given the n available objects searched (see Table 3).

Figure 5 presents the line fluxes and their significance for the UV line detections together with a sample of literature values (small dots; see Appendix C for details). For the emission line doublets the total flux from combining both components is shown. The uncertainties are determined by the variance of the template crossmatch defined in Eq. (A.5). The large number of upper limits provided by the FELIS template matching are not shown to prevent cluttering the figure, but are available in the public catalog published with this paper. As a crosscheck, we compared the estimated He II and C IIIemission line fluxes from FELIS to the PLATEFIT (Tremonti et al. 2004; Brinchmann et al. 2004) flux estimates for the sources in the UDF mosaic also studied by Nanayakkara et al. (2019) and Maseda et al. (2017) and find good agreement between the fully independent detections and flux measures.

Fig. 5. Distribution of rest-frame UV emission line fluxes (x-axis) detected at a FELIS S/N (y-axis) above 3σ (horizontal dashed line) in the MUSE-Wide (filled circles), the UDF mosaic (filled diamonds), and the UDF10 (filled x’s) data sets studied in this paper. These detections are compared to the collection of detections from the literature (small dots) described in Appendix C. Here S/N represents the emission line S/N quoted in the original references. All points are color coded according to their redshift. The solid histograms show the subset of the objects from this work. The dotted histograms show the distribution of measurements when including the literature data.

Figure 6 presents comparisons between the UV emission line fluxes and the C III flux for objects with both C III and additional UV lines detected. It is clear that several of the emission line fluxes correlate in an approximately linear fashion. To quantify the strength and parametrize these correlations we estimate the Pearson correlation coefficient (r_P) and the Spearman’s rank correlation coefficient (r_S). The r_P tests the strength of a linear relation, whereas r_S only requires a monotonic correlation to provide a larger value. The coefficients have values between −1 and 1 and as a rule of thumb absolute values of the correlation coefficients |r|< 0.3, 0.3 < |r|< 0.5, 0.5 < |r|< 0.7, and |r|> 0.7 can be considered very weak (nonexistent), weak, moderate, and strong correlations, respectively. Standard linear regression usually does not account for uncertainties in both data sets. To obtain a parametric representation of the correlations, we therefore estimate the best-fit linear relations between the flux estimates via orthogonal distance regression (ODR) using Scipy’s build-in version of this. The correlation coefficients and resulting best fits for the emission line fluxes are presented as correlations 1−4 in Table 4. The correlations between C III and He II, O III, and Si III are tighter (higher values of r_P) than the correlation with the resonant emission from C IV. In all cases, the upper 3σ limits resulting from nondetections of UV emission lines by the FELIS template matches follow the presented correlations. Hence, if the redshift evolution is insignificant (no indication of a redshift dependence is seen for the data in the redshift range probed by the MUSE samples) these correlations can be used to predict likely emission line fluxes for the fainter UV emission features based on C III. This will be useful for targeting these fainter lines in the EoR where only the brightest UV lines (Lyα and in some cases C III and C IV) have been recovered with current facilities.

Fig. 6. Similar to Fig. 5, but showing correlations between the > 3σ C III detections and the other secondary UV emission line fluxes from MUSE-Wide, the UDF mosaic, and UDF10 shown as large symbols. The literature measurements described in Appendix C are shown by the small dots. The dashed lines show the one-to-one relations to guide the eye.

If we instead consider the individual components of the UV emission line doublets, we find that the flux ratios F(N V₁)/F(N V₂) and F(O III]₁)/F(O III]₂) are both generally < 1, whereas F([Si III])/F(Si III]) and F([C III])/F(C III]) are > 1. The mean, median and standard deviation of these four flux ratios are (0.69,0.65,0.24), (0.49,0.44,0.25), (1.65,1.60,0.73) and (1.35,1.37,0.45). The data are shown in Fig. B.1. The scatter for the F(C IV₁)/F(C IV₂) correlation is larger (lower r_P) than for the nebular emission lines and the values of this flux ratio have mean, median and standard deviation (1.34,1.30,0.83). This can be attributed to the line profiles of the C IV emission which are poorly captured by the Gaussian FELIS line templates. C IV is a resonant line that is scattered by high-ionization gas and is produced in the winds of young massive O and B stars as well as in the ISM (e.g.; Leitherer et al. 2011; Du et al. 2016; Berg et al. 2019a; Feltre et al. 2020). The wind features can produce prominent P Cygni profiles which, superimposed on the ISM emission, create line-profiles which are poorly approximated by a Gaussian. Hence, when such a combination of emission processes is present the flux measures (and estimates of velocity offset with respect to systemic; see Sect. 8) are more uncertain for the C IV template matches and are therefore expected to scatter more. Full modeling of the C IV line profiles is beyond the scope of this work, but provides valuable information about the underlying emission mechanisms, the metallicity and the initial mass function Leitherer et al. (2011).

Even though, an attempt to fully model the stellar+nebular emission lines like, He II, C IV, and N V is beyond scope of this work, broader lines measured from the FELIS template matches, could be indicative of a predominantly stellar contribution to the emission (for non-AGN), as stellar emission is generally broader than the nebular emission (e.g., Shirazi & Brinchmann 2012; Crowther et al. 2006; Crowther 2007; Nanayakkara et al. 2019). The FWHM distributions from the FELIS matches for these lines span the range between roughly 50 km s⁻¹ and 550 km s⁻¹ with the majority of detections having FWHM ≲ 300 km s⁻¹ for all three lines. Thus, the FWHM of He II, C IV, and N V are generally modest, but could indicate contribution from stellar emission in some systems. For details on the nature of the HeII emission and the emitter properties we refer to Nanayakkara et al. (2019), who as mentioned studied the majority of the HeII emitters presented here in detail.

By comparing and fitting literature compilations of the oscillator strength Martin & Wiese (1976), Feibelman (1983) and Morton (1991) predict the theoretical C IV doublet ratio to be F(C IV₁)/F(C IV₂) = 2 as listed in Table 1. However, doublet ratios closer to one have also been found (e.g., Feibelman 1983; Christensen et al. 2012; Stark et al. 2014). The F(N V₁)/F(N V₂) emissivity ratio is also expected to be two (Martin & Wiese 1976; Torres-Peimbert & Pena 1984; Morton 1991). This is in disagreement with the N V flux ratios measured by FELIS for the seven MUSE sources with potential N V detections, where we generally see that F(N V₂) > F(N V₁). As N V is also a resonant line arising from stellar winds and in the ISM (Leitherer et al. 2011), profiles deviating from Gaussian could be part of the explanation for this discrepancy. The best-fit linear correlations and correlation coefficients for these two emission doublets are presented as correlations number five and six in Table 4.

As listed in Table 1Morton (1991) estimates the O III doublet ratio to be F(O III]₁)/F(O III]₂) ≈ 0.7 which is in agreement with recent findings where O III]₂ tends to be strongest (e.g., Mainali et al. 2017; Vanzella et al. 2016). This is also in agreement with what we find here. Performing a linear fit to the measured O III doublet component fluxes from our sample and the literature and estimating the correlation coefficient we find the relation presented as correlation seven in Table 4. Assuming an electron temperature of 10⁴ K the theoretically expected Si III doublet ratio is F([Si III])/F(Si III]) ≲ 1.7, similar to the expected C III doublet ratio F([C III])/F(C III]) ≲ 1.6 (Osterbrock & Ferland 2006). As we show in Sect. 9, these two ratios depend somewhat on the assumed electron temperature and probe the electron density of the emitting gas. As mentioned, we find that both of these flux ratios are > 1 for the vast majority of sources. The best-fit ODR linear empirical relations and correlation coefficients for these ratios are presented as correlations number eight and nine in Table 4 and are in agreement with the theoretical expectations.

6. Rest-frame EW estimates

In addition to the line flux measurements of the UV emission line detections from FELIS, we compute the rest-frame equivalent width (EW₀) using the expression

(2)

The UV emission line flux (F_{UV emission,  rest-frame}) is provided directly from the FELIS template matching of the emission line templates (α in Eq. (A.4)). The observed continuum flux density (f_{continuum,  observed}) is estimated from available photometric catalogs in the following way: For each source we identify the HST broad-band nearest to the location of the detected UV emission line free of any emission line contamination. We then assume a power law continuum of the form f(λ)∝λ^β with a fixed spectral slope of β = −1.97, which is the median spectral slope for the LAEs studied here and by Kerutt et al. (2021). If available, we use the photometry by Kerutt et al. (2021, only estimated for LAEs) described in Sect. 7 to predict the continuum flux density at the location of the relevant UV emission. Otherwise we use the Rafelski et al. (2015) or the Skelton et al. (2014) HST broad-band fluxes to estimate the continuum level. Assuming β = −1.97 for the non-LAE (z < 2.9) is also a good approximation of the average spectral slope based on the available HST photometry of these objects in the MUSE wavelength range. Figure 7 shows the HST broad-band magnitudes used to infer the continuum flux density for each of the detected UV emission lines for sources with estimated EW₀.

Fig. 7. AB magnitude distributions of the measured HST magnitudes used to infer the underlying continuum flux density when estimating the EW0 of the UV emission lines targeted in this study.

In Fig. 8 we show the EW₀ values resulting from a combination of these continuum flux density estimates and the emission line fluxes from the FELIS template matches. We note that the apparent “correlations” seen in these panels are driven by the fact that the continuum magnitude shown on the x-axes is included in the definition of the EW₀ estimates plotted on the y-axes. Taking out the dependence on the continuum magnitude of the EW₀, we see that the line fluxes as a function of continuum AB magnitude are fairly flat (with a large scatter), with a tendency for fainter lines to correspond to fainter objects. The full sample of EW₀ estimates (and upper or lower limits; not shown in Fig. 8) are available in the value added catalog described in Appendix B.

Fig. 8. Rest-frame equivalent width (EW0) of the UV emission lines detected by FELIS as a function of the continuum magnitude used to estimate the flux density at the emission line location. As in Fig. 5, filled circles, diamonds and x’s correspond to objects from the MUSE-Wide, the UDF mosaic and the UDF10 samples, respectively. Each point is color coded according to each object’s redshift.

As can be seen in Fig. 8, a few sources show high EW estimates that are rarely seen in the literature, For example, we see objects predicted to have EW(C III) = 30−300 Å and EW(C IV) above 40 Å, and a single source with EW(He II) ≈ 100 Å. If any of the high-EW measurements discussed below are confirmed, these systems will provide interesting extremes for further studies. The high-EW(C III) systems are 603092083, 301006546, and 721870849. The potential detection from object 603092083 coincides with a sky feature and its residuals, which have likely affected the measured line fluxes. Object 301006546 has a potential detection of both C III and C IV. Both of these detections are however low-significance detections and the large EW estimate questions whether these detections are real. Finally, the prominent C III emission from object 721870849 is present in a spectrum with no apparent continuum presumably arising from a low-luminosity system resulting in the extreme EW at z = 2.3796. The photoionization models discussed in Sect. 11 provide non-AGN solutions that are capable of reproducing the observed flux ratios of this system, though the parameter space is fairly limited. The object with a potential He II EW around 100 Å (object 601071350) also has EW(C IV) above 40 Å (which is however affected by sky residuals). As pointed out by Nanayakkara et al. (2019), who also studied this object (ID = UDF 3621), a nearby AGN and an LAE can affect the measured HeII emission in this system. Our flux measurements agree with those from Nanayakkara et al. (2019), whereas our EW estimate deviates, as they find EW(He II) = 6.8 Å for this object. The three remaining objects with EW(C IV) > 40 Å not discussed above are 102014087, 210012237, and 721480767. The estimated EW(C IV) of 58 ± 8 Å, 45 ± 10 Å, and 46 ± 6 Å from the spectra of these objects are generally of lower quality with S/N(FELIS) of 3.6, 3.4, and 4.3, which could indicate that the detected lines could be spurious or their strengths less certain.

Considering correlations between the EWs of the various UV emission lines, similar to what is presented for the emission line fluxes in Fig. 6, reveals that EW₀(Si III), EW₀(He II), and EW₀(O III) generally follow the strength of EW₀(C III), though offset from the one-to-one relation (as shown in Fig. B.2). The C IV resonant line does not correlate with EW₀(C III) to the same degree as the nonresonant UV lines. As there are only seven objects with both C III and N V detected we cannot draw any firm conclusions but no clear correlation appears to be forming. The linear ODR fits to the logarithmic distributions of the EW₀ estimates from the current study and the literature are presented together with the correlations coefficients in Table 4 as correlations number 10−14. Again, the limits for the nondetections are in agreement with the presented correlations. Hence, similar to the correlations for the emission line fluxes found above, EW₀(Si III), EW₀(He II), EW₀(O III), and EW₀(C III) correlate with each other and can be used as predictive tools for estimating EW₀ of emission lines. Based on photo-ionization models, a correlation between the EW₀ of C III and O III is expected as shown in Fig. 9 of Jaskot & Ravindranath (2016).

In line with these findings, multiple studies including Maseda et al. (2017), Mainali et al. (2020), and Tang et al. (2021a) have reported that high-EW C III emitters generally have high EWs of [O III]λ5007 + Hβ as also anticipated by the Jaskot & Ravindranath (2016) theoretical models. Jaskot & Ravindranath (2016) stress that both EW₀([O III]λ5007 + Hβ) and EW₀(C III) reach their largest values for young, high ionization parameter models, but confirm that C III emission is more sensitive to metallicity due to its temperature dependence. Tang et al. (2021a) show that EW₀(C III) depends strongly on the metallicity of the emitting system and that you generally need low metallicity to obtain large EW₀(C III). They find that EW₀(C III) increases by a factor three when the metallicity changes from 0.3 to 0.1 of solar metallicity, whereas [O III]λ5007 + Hβ varies little with metallicity in their models. Hence, Tang et al. (2021a) argue that to obtain EW₀(C III) > 10 Å requires an ionization parameter log₁₀(U) above −2.5 and a metallicity of 0.2 solar or lower. The fact that we see several such systems at redshift three (lower right panel in Fig. 8) indicates that even at redshift 3, the large sample of MUSE sources includes low-metallicity objects that are typical for objects at redshifts approaching the EoR. Furthermore, it was pointed out by Chevallard et al. (2018) that systems with large EW₀([O III] λ4959 + [O III] λ5007) and hence large EW₀(C III) tend to have a high production efficiency of H-ionizing photons, ξ_ion.

7. The UV emission of the Lyα-emitter sample

In this study we are focusing on rest-frame UV emission red-wards of the strongest UV emission line, Lyα, partially motivated by the challenges of observing this line at high redshift (z ≳ 6) where the significantly neutral CGM and IGM absorbs the Lyα photons escaping the galaxy along the line of sight (e.g., Dijkstra et al. 2011; Laursen et al. 2011, 2019; Dijkstra 2017). Nevertheless, the Lyα line itself has improved our understanding of star-forming galaxies in the (early) Universe. In particular, the asymmetric Lyα line profile has enabled redshift confirmations of large samples of sources at both 2 < z < 6 (e.g., Steidel et al. 2014; Le Fevre et al. 2015; Herenz et al. 2017; Inami et al. 2017; Urrutia et al. 2019) and high redshift at z > 6 (e.g., Finkelstein et al. 2013; Oesch et al. 2015; Schmidt et al. 2016; Tilvi et al. 2016; Huang et al. 2016; Pentericci et al. 2018; Fuller et al. 2020). The resonant scattering of the photons and the resulting (occasional) multipeaked emission has been shown to relate closely to the column density and dynamics of the neutral hydrogen in the ISM and the CGM (Verhamme et al. 2015; Gazagnes et al. 2018, 2020). The fraction of galaxies with confirmed Lyα emission has been used to probe the evolution (or lack thereof) of the fraction of LAEs among Lyman-break galaxies from low redshift to the EoR (e.g., Treu et al. 2013; Pentericci et al. 2014; Tilvi et al. 2014; de Barros et al. 2017; Caruana et al. 2018; Kusakabe et al. 2020). Together with the observed velocity offset of the Lyα line resulting from resonant scattering (Schenker et al. 2013; Erb et al. 2014; Hashimoto et al. 2015; Stark et al. 2017; Verhamme et al. 2018), this has probed the amount of neutral gas in the IGM and has constrained the neutral fraction of the Universe during the EoR (Ouchi et al. 2010; Greig et al. 2017; Mason et al. 2018a,b; Banados et al. 2018; Hoag et al. 2019). Furthermore, comparisons between Lyα and Hα or UV emission line strengths have been used to study the production efficiency and escape of ionizing photons from LAEs (Nakajima et al. 2016; Matthee et al. 2017; Harikane et al. 2018; Lam et al. 2019; Maseda et al. 2020). It is therefore of interest to relate and compare the measured rest-frame UV emission lines red-wards of Lyα studied here, with the characteristics of the Lyα line itself and the properties of the LAEs in our sample.

For this comparison we rely on HST broad-band magnitudes, Lyα line fluxes, Lyα EW₀ estimates, spectral UV slopes β, continuum magnitudes, and Lyα FWHM from the catalog that will be presented by Kerutt et al. (2021). This LAE study is based on the same data and source identification (see Sect. 3) as the current study, but focuses on properties of the emanating Lyα emission. The Lyα emission fluxes correspond to the measured flux within 3D apertures of three Kron (1980) radii as measured by LSDCat when detecting sources in the MUSE data cubes. We use these Lyα fluxes as opposed to obtaining them directly from the TDOSE spectra, as the TDOSE extractions are based on the assumption that the morphological extent of the line emission follows the continuum morphology of the modeled HST images (Sect. 4 and Schmidt et al. 2019). However, Lyα emission is known to be extended beyond the continuum (Steidel et al. 2011; Momose et al. 2014; Wisotzki et al. 2016, 2018; Leclercq et al. 2017, 2020) and fluxes based on the TDOSE spectra would therefore be biased. The LSDCat Kron radii fluxes therefore better represents the actual Lyα flux emitted by the LAEs. As for the secondary UV emission lines, the EW₀(Lyα) values were calculated by comparing the fluxes to the continuum flux densities estimated from a continuum represented by a power law f_λ ∝ λ^β. To obtain the spectral slope, Kerutt et al. (2021) first determined the magnitudes from available ancillary broad-band HST photometry, by fitting the rest-frame UV morphology for each of the LAEs using GALFIT (Peng et al. 2010, 2002). This provided morphological parameters for all LAEs including a measure of their effective radius. The estimated absolute UV magnitude at 1500 Å is also based on these GALFIT models. The spectral slope β was then obtained from fitting the continuum power law to the GALFIT-based HST magnitudes. To avoid large scatter in the EW₀ measurements presented and analyzed by Kerutt et al. (2021), the EWs are based on the median β for the full LAE sample of β = −1.97 similar to what was done for the secondary UV emission lines presented here. The FWHM of the Lyα emission was measured for each source by fitting a skewed Gaussian profile (Eq. (2) by Shibuya et al. 2014) to the Lyα line profiles in 1D spectral extractions weighted by the MUSE PSF to maximize S/N. These fits also provide a Lyα redshift which is more precise than the lead line redshits provided by the LSDCat source identification. We therefore use these redshifts for the analysis of the Lyα velocity offsets described in Sect. 8. Finally, Kerutt et al. (2021) provide estimates of the systemic redshifts based on the FWHM and peak separation between any double-peaked LAEs (identified through visual inspection of the 1D spectra) in the sample based on the empirical relations presented by Verhamme et al. (2018). We note that a handful of the z > 2.9 objects studied here are not included in the Kerutt et al. (2021) catalog, as their selection was based only on non-AGN objects with leading Lyα emission based on the LSDCat selections (Sect. 3). Hence, for z > 2.9, the objects with IDs 121033078, 601381485, 720470421, 722551008, 722731033, and 723311101 are not included in the LAE parameter comparisons in the following. For further details and for the full value-added catalog of Lyα-related quantities we refer to Kerutt et al. (2021).

In Fig. 9 we compare EW₀(Lyα) from the Kerutt et al. (2021) catalog to the estimated UV emission line fluxes detected with FELIS in the 1997 LAEs (z > 2.9). We see no strong correlations in general. It does however appear that the objects with larger EW₀(Lyα) mostly have C III or O III UV emission lines detected at larger fluxes, that is there are no EW₀(Lyα) > 200 Å objects with C III or O III line detections below ≈2 × 10⁻¹⁸ erg s⁻¹ cm⁻². The few Si III detections from the current study and available from the literature only arise in sources with EW₀(Lyα)≲100 Å.

Fig. 9. Detected emission line fluxes as a function of rest-frame equivalent width of Lyα. The measurements presented here from MUSE-Wide, the UDF mosaic and the UDF10 are shown as large symbols. Measurements from the literature are shown by small dots (see Appendix C). All symbols are color coded according to source redshifts and are shown with 1σ error bars. The solid histograms show the subset of the objects from this work, whereas the dashed histograms include the estimates from the literature collection.

Figure 10 presents the current census of estimated UV emission line EW₀ from our study (large symbols) and the literature (small dots, see Appendix C) for LAEs as a function of EW₀(Lyα). For reference, we show lines of fixed ratios between the UV emission line EW₀ and EW₀(Lyα) estimates.

Fig. 10. UV emission line EW0 estimates of LAEs as a function of EW0(Lyα). The measurements from the current study are shown as large symbols (circles from MUSE-Wide, diamonds from the UDF mosaic, and x’s from UDF10). Measurements from the literature are shown by small dots (see Appendix C). All symbols are color coded according to source redshifts and are shown with 1σ error bars. The diagonal curves correspond to the relations 10:1, 3:1, 1:1, 1:3, and 1:10 for reference. Solid histograms show objects presented in this study, whereas dotted histograms include the literature samples.

Similar to the correlations of the EW₀ estimates for the UV emission lines probing systemic redshift reported in Sect. 6, we see that several of the UV lines have strengths that can be related to the Lyα EW. Except for the resonant C IV emission, the N V and O III emission, which all have correlation coefficients of roughly 0.5 or below, the UV lines show fair correlations with EW₀(Lyα). For completeness, in Table 4 we list all six ODR fits to the data shown in Fig. 10. The majority of the EW₀(C III) measurements appear to align just below the dashed line marking EW₀(C III) of one-third the strength of the Lyα lines for EW₀(Lyα) > 10 Å. If we focus on the subset of objects with EW₀(C III) < EW₀(Lyα) the distribution of the ratio between the two EWs has a mean value of 0.22 ± 0.18. For the full sample the mean ratio is 0.36 ± 1.46. Similar trends, in some cases with equally large amounts of scatter, have been presented in various other studies (e.g., Shapley et al. 2003; Stark et al. 2014, 2015a; Rigby et al. 2015; Du et al. 2018; Le Fevre et al. 2019; Marchi et al. 2019). As noted in Sect. 6, the strength of C III depends on the metallicity of the galaxy but is also affected by the hardness of the ionizing radiation (Erb et al. 2010; Senchyna et al. 2017). This has been shown to lead to expected correlations between EW₀(C III) and EW₀(Lyα) in photoionization models (Jaskot & Ravindranath 2016; Harikane et al. 2018). Hence, despite that the strength and significance of the observed correlations between EW₀(C III) and EW₀(Lyα) are still being debated in the literature, the correlation presented here (correlation 20 in Table 4) appears to be in agreement with both previous observational results and theoretical expectations.

In addition to assessing potential correlations between the EW₀ and emission line fluxes derived in this work, we also considered correlations with the LAE spectral slope β, Lyα line peak separation, the apparent and absolute magnitudes at 1500 Å (m_UV, 1500 and M_UV, 1500), LAE effective radius, (systemic) redshift, Lyα flux, Lyα luminosity and FWHM(Lyα), all from the Kerutt et al. (2021) LAE property catalog. Overall, we did not find prominent trends between these parameters and line fluxes and EWs of the secondary UV emission line with the exception of tentative correlations with the Lyα luminosity and the spectral slope β.

We find correlations between the Lyα luminosity and the flux of the detected UV emission lines C III, He II, N V, and O III with 0.62 < r_P < 0.83. When accounting for the continuum level at the location of the emission features via the EW₀ estimates, the tentative correlations persist (−0.72 < r_P < −0.64) except for the O III measurements. For the objects with detected C III a larger EW₀(C III) seems to imply a lower Lyα luminosity when 41 < log(L_Lyα) < 43. The opposite appears to be true for objects with He II and N VḢowever, this is based on only a handful of objects in each case, and is therefore prone to scatter in the measurement and their uncertainties. For all lines the UV emission line flux increases with the Lyα luminosity. Figures B.3 and B.4 present the Lyα luminosity, the UV emission line fluxes and EW₀ estimates.

The F(O III) and EW₀(O III) of the LAEs where O III is detected correlates with the spectral slope β with r_P = 0.97 and −0.82, respectively. This indicates that O III-emitters with bluer (β ≈ −2.1) spectral slopes have larger EW₀(O III) but smaller F(O III) than objects with redder slopes (β ≈ −1.6). A possible explanation for this trend is that bluer galaxies are likely younger and therefore have harder ionizing spectra leading to stronger UV lines. This agrees with the trends seen with UV emission lines in the composite MUSE UDF mosaic spectra presented by Feltre et al. (2020). For the other UV lines the correlation coefficients r_P and r_P are all below 0.46 and therefore provide no clear indication of potential trends. For this comparison we used the individually measured β values from Kerutt et al. (2021). The potential correlations between O III and the spectral slope β are shown in Fig. B.5.

8. Emission line velocity offsets

As noted, one of the key diagnostics obtainable from the Lyα line is the velocity offset of the resonant Lyα emission, Δv_Lyα, from the systemic redshift. In more general terms the offset of the Lyα emission, or any other line, can be measured with respect to any reference (emission line) wavelength whether it is at systemic or not. For instance, the relative velocity offset between the resonant C IV emission and the Lyα line could potentially reveal similarities or differences between the neutral hydrogen that scatters Lyα photons and the ionized gas that scatters the energetic C IV photons. Or the offset of C IV with respect to systemic can be used as a tracer of the ionized gas’ velocity structure similar to what has been done for Lyα and neutral hydrogen. Following Erb et al. (2014) and Shibuya et al. (2014), we therefore define the general emission line velocity offsets as

(3)

where c is the speed of light in km s⁻¹, “line” refers to the emission line for which the velocity offset is measured, and “reference” refers to the reference feature used to measure the offset. Hence, the reference redshift, z_reference, is either the systemic redshift, or the redshift of the observed location of another line.

In Fig. 11 we present the UV emission line velocity offsets from the FELIS detections between the LSDCat lead emission line and the respective UV emission lines. For z ≳ 2.9, which is marked by the vertical dashed lines, this effectively means that the lead line is Lyα and hence this part of each panel shows Δv_Lyα. For the LAEs we replace the LSDCat redshift with the Lyα redshift from the skewed Gaussian Lyα profile fits when determining the velocity offsets. When estimating the velocity offsets, we did do not impose any cut on S/N(FELIS). As illustrated by the color coding in Fig. 11 the S/N(FELIS) values for the detections used to estimate velocity offsets mostly fall around S/N ≈ 5 (green). When visually inspecting the potential line detections from FELIS, we made an explicit cut to only consider emission features with velocity offsets below 1000 km s⁻¹ with respect to the primary line of the object. As AGN and quasars are known to emit rest-frame UV lines at velocity offsets up to several thousands km s⁻¹ (e.g., Onoue et al. 2020), our limitation to only consider candidate detections with Δv < 1000 km s⁻¹ potentially limits our ability to detect and recover emission from faint AGN in our sample. However, this limitation does not prevent us from recovering several known AGN in the targeted fields with velocity offsets below 1000 km s⁻¹ (see, for example, Fig. E.1). In cases where one component of a detected emission line doublet is not detected the solution with the highest S/N(FELIS) from the simultaneous doublet template match (see Sect. 5) is used to obtain the systemic redshift. Of the ancillary spectroscopic redshift collected in support of the photometric redshifts presented by Skelton et al. (2014), none are of the sources presented with systemic redshifts estimates from MUSE presented here.

Fig. 11. Velocity offsets for the individual FELIS template matches of the detected UV emission lines with respect to the object redshift based on the “lead line” of the MUSE LSDCat redshift determinations. For LAEs the Lyα redshift based on the line profile fits by Kerutt et al. (2021). From top left to bottom right the velocity offset for C IV, He II, O III, Si III, and C III with respect to the lead line are shown. For objects with z ≳ 2.9 (vertical dashed line) the lead line of essentially all objects is Lyα. The color coding shows the significance of the FELIS detection in terms of the signal to noise ratio, S/N(FELIS). All Δv estimates were assigned an uncertainty of ±100 km s−1 according to the precision uncertainty coming from the MUSE spectral resolution according to the discussion in Appendix A.3.

Figure 11 shows that for the tracers of the systemic redshift with good statistics (He II, O III, and C III) the Lyα velocity offsets are 250−500 km s⁻¹ for the majority of the objects. For non-LAEs on the other hand, most objects are consistent with no velocity offsets within 3σ. These Δv_Lyα estimates agree well with what has been found in the literature at both low and high redshift (e.g., Rakic et al. 2011; Erb et al. 2014; Stark et al. 2017; Mainali et al. 2017; Cassata et al. 2020).

There are two exceptions to this trend among these LAEs. For object 604992563 the velocity offset between Lyα and C III is estimated to be −145 km s⁻¹. This either questions the potential detection of the C III doublet (the FELIS S/N for this detection is estimated to be 4.0) or could indicate infalling gas offsetting the Lyα blue-wards of systemic (e.g., Verhamme et al. 2006; Dijkstra et al. 2006; Mitchell et al. 2021). However, within 2σ the estimated offset is still consistent with zero or a small positive offset of the Lyα line. We show the spectrum of object 604992563 in Fig. E.2. The second LAE showing significant apparent blueshift of the Lyα emission is object 219009247 which has a potential He II emission detected −634 km s⁻¹ offset from the Lyα redshift. Upon further inspection, it turns out that the emission detected as He II for object 219009247 is actually [O III] λ5007 flux from a chance superposition (projected distance ) of a foreground line emitter at z = 0.3731 contaminating the spectrum of object 219009247. Hence, this detection presents a false positive detection of He II emission from a LAE in our sample.

Empirical correlations between Δv_Lyα and EW₀(Lyα) (see Tang et al. 2021b, for a recent collection) have been parametrized through empirical relations by, for example, Nakajima et al. (2018) and Adelberger et al. (2003) for LAEs at z ≈ 3. In Fig. 12 we show the Δv_Lyα estimates for our sample based on the C III detections together with these empirical relations. The MUSE data appear to roughly follow the Nakajima et al. (2018, dot-dashed) relation though with large scatter. This amount of scatter is however comparable to the scatter of the collection of data that Nakajima et al. (2018) based their relation on.

Fig. 12. Velocity offsets of Lyα determined with respect to the C III detections as a function of the Lyα rest-frame equivalent width estimated for a fixed UV slope of β = −1.97. The gray dot-dashed and dotted lines show the empirical relationship between ΔvLyα and equivalent width for z ≈ 3 objects presented by Nakajima et al. (2018) and Adelberger et al. (2003), respectively. The individual points are color coded according to redshift.

Previous studies, among those Erb et al. (2014), have noted a correlation between absolute UV magnitude and Δv_Lyα. In Fig. 13 we show the estimated UV magnitudes for the Lyα velocity offsets of the C III emitters in our sample (large symbols) together with a collection of M_UV and Δv_Lyα measurements from the literature shown as small dots (see figure caption for details). All points have been color-coded according to the object redshift and they indicate that the highest redshift objects (blue and purple points) on average have smaller Lyα velocity offsets and are brighter than objects at redshifts 2−3 (green points). The latter could of course be affected by the Malmquist (1920, 1922) bias, as it is generally harder to observe intrinsically faint objects at higher redshifts. For comparison, we show a set of colored dashed curves predicting the correlation between the median Δv_Lyα and M_UV from Mason et al. (2018b). These curves correspond to the black curve in their Fig. 2 at different redshifts (according to the color coding) which relates the velocity offsets to the galaxy halo mass (their Eq. (2)). The galaxy halo masses are translated into absolute UV magnitudes through a simple abundance matching model. Based on the relation between M_UV and halo mass of the emitting systems by Mason et al. (2015); Mason et al. (2018b) suggest that Δv_Lyα probes the halo mass of the host galaxy (though with significant scatter), as the amount of neutral hydrogen scattering Lyα is closely related to the galaxy mass. Despite the significant scatter in the points plotted in Fig. 13 they appear to roughly follow the curves from Mason et al. (2018b). This could indicate that the halo mass, that is the reservoir of available gas to scatter the Lyα photons, is indeed a more important quantity for determining the Lyα velocity offsets than, for instance, star-formation rate (SFR), that is outflows from star formation. If the opposite was true, at fixed M_UV the high-redshift galaxies should show larger Δv_Lyα than the lower redshift galaxies as galaxies generally have lower mass for a fixed SFR at higher redshifts. The lowest redshift targets from the literature (z ≲ 1.5), that is the yellow points, seem to be an exception to this apparent trend, potentially indicating that for these systems the hypothesis that the halo mass drives the size of the velocity offsets might not be true. Here star formation processes might be of higher importance.

Fig. 13. Velocity offsets of Lyα as a function of the absolute UV magnitude for the MUSE LAEs with C III detections studied here shown by the large symbols and a collection of measurements at redshifts from z ≈ 2 to z = 7.73 collected from Kulas et al. (2012), Schenker et al. (2013), Erb et al. (2014), Stark et al. (2015a, 2017), Henry et al. (2015), Willott et al. (2015), Pentericci et al. (2016), Inoue et al. (2016), Bradač et al. (2017), Mainali et al. (2017), Verhamme et al. (2017), Sobral et al. (2018), Hutchison et al. (2019), Cassata et al. (2020), and Matthee et al. (2020a,b), shown by the small dots. To prevent cluttering the figure the estimated uncertainties are not shown. The median errors on the magnitudes and ΔvLyα are 0.2 mag and 50 km s−1. The dashed curves present the predicted median ΔvLyα as a function of UV magnitude for a range of redshifts (according to the color coding) from Mason et al. (2018b, corresponding to the black curve in their Fig. 2 from their Eq. (2)). This is based on an abundance matching model relating the absolute UV magnitude to the galaxy halo mass (Mason et al. 2015). Points and curves are color-coded according to object redshift and the color gradient highlights the importance of halo mass, that is available gas reservoir, in producing and determining ΔvLyα. The solid histograms show the distribution of objects from our study, whereas the dotted histograms include the sample of measurements from the literature.

Marchi et al. (2019) highlight the importance of ISM outflows in producing velocity offsets of Lyα. However, in line with the above, they argue that the key factor controlling the strength of the outflows, and hence the Lyα velocity offsets, is the H I column density. Only in systems with low column density (≈10¹⁹ cm⁻²) can the ISM outflows produce velocity offsets, which are not expected to be larger than ≈300 km s⁻¹. This appears to be in agreement with the literature z ≲ 1.5 objects which all have Δv_Lyα ≲ 300 km s⁻¹. According to Marchi et al. (2019), to produce larger Δv_Lyα it appears that systems with larger gas reservoirs (halo masses) with higher H I column densities are needed. Based on samples of star forming galaxies at z ≈ 0.2 Shibuya et al. (2014) and Henry et al. (2015) arrive at similar conclusions, which are also supported by radiative transfer models (e.g., Chonis et al. 2013; Verhamme et al. 2015).

Muzahid et al. (2020) estimated the average velocity offset of Lyα for a sample of LAEs from the MUSE Quasar-field Blind Emitters Survey (MUSEQuBES; PI. Schaye; Chen et al. 2019) with respect to circumgalactic absorption lines, and found an average of Δv_Lyα = 171 ± 8 km s⁻¹, which appears to be somewhat lower than what we find here. Their sample generally contains objects with smaller L(Lyα) than the UV line emitters providing Δv_Lyα here, which could explain part of this difference as brighter Lyα implies larger halo mass and brighter M_UV (e.g., Kerutt et al. 2021; Khostovan et al. 2019) which again implies larger Lyα velocity offsets at fixed redshift (Mason et al. 2018b).

Based on another sample of Lyα velocity offsets, Verhamme et al. (2018) present means of predicting the systemic redshifts based on the FWHM of the Lyα line or the separation of the red and blue components of double-peaked Lyα emission profiles. Figure 14 presents a comparison of the Lyα velocity offsets for the C III emitters among the LAEs in our sample and the predicted Lyα velocity offset based on the Verhamme et al. (2018) relations. Earlier measurements on five sources from our main sample were included in the work by Verhamme et al. (2018). If a Lyα profile was classified as double peaked, the correlation between velocity offset and peak separation was used. Otherwise the correlation with Lyα FWHM was used to predict the velocity offset. Figure 14 shows that there is good agreement between the measured and predicted Lyα velocity offsets. Similarly, the Lyα velocity offsets based on the He II and O III detections also match the empirical predictions presented by Verhamme et al. (2018).

Fig. 14. Comparison of the velocity offset of Lyα with respect to the C III doublet detections and the Lyα velocity offset estimated based on the peak separation for double peaked LAEs and the Lyα FWHM for single peaked LAEs following the empirical relations presented by Verhamme et al. (2018). The dashed line indicates the one-to-one relation.

If we consider the two individual predictions from the original study shown in Fig. 15 relating Δv_Lyα to the Lyα peak separation (top) and FWHM (bottom), we see that our measurements (large symbols) are in good agreement with the empirical relations from Verhamme et al. (2018):

Fig. 15. Comparison of the velocity offset of Lyα and half the Lyα peak separation for objects with Lyα line profiles containing both a red peak and a blue bump (top) and the Lyα FWHM (bottom). These are the quantities used to predict ΔvLyα (Verhamme et al. 2018 approx.) on the x-axis in Fig. 14. The gray bands show the empirical relations presented by Verhamme et al. (2018). The small dots present a collection of literature measurements from Kulas et al. (2012), Leitet et al. (2013), Oesch et al. (2015), Hashimoto et al. (2015), Henry et al. (2015), Guaita et al. (2016), Izotov et al. (2016a,b), Stark et al. (2017), Verhamme et al. (2017, 2018), Yang et al. (2017), Matthee et al. (2020b), Muzahid et al. (2020), and Matthee et al. (2021). All points are color coded according to object redshift. The solid histograms show the measurements presented in this study, whereas the dotted histograms include the measurements from the literature.

These relations are shown as the gray bands in Fig. 15 and the correlation coefficients for each of them are estimated after including the collection of measurements from the literature shown as the small symbols in the figure. Our measurements appear to be offset slightly high with respect to the Δv_Lyα − FWHM(Lyα) relation, even though we also see that they present a scatter similar to what is found in the literature. As an estimate of the scatter in the relations, we determine the perpendicular Euclidian distance between all points and the proposed linear relations. The average distance from the peak separation (FWHM) relation is 70 km s⁻¹ (68 km s⁻¹) with a standard deviation of the distribution of distances of 77 km s⁻¹ (67 km s⁻¹). From the data used to define the empirical relations Verhamme et al. (2018) quote a similar but tighter scatter of 53 ± 9 km s⁻¹ and 72 ± 12 km s⁻¹ for the two relations and determine an uncertainty on the systemic redshift correction of ≲ ± 100 km s⁻¹ which is in agreement with the above estimates.

In addition to the correlations described above, we also checked for trends between Δv_Lyα and the remaining LAE properties described in Sect. 7, including β, LAE effective radius, Lyα flux, and L(Lyα). We do not find any trends between these parameters and the Lyα velocity offset with respect to the UV emission probes of systemic redshift (|r_P|< 0.42 and |r_S|< 0.28 for all parameters). The only exception to this is the potential correlation (shown in Fig. B.6) between Δv_Lyα and the effective radius indicating that larger galaxies (R_e ≈ 3.5 kpc) have offsets of roughy 500 km s⁻¹ as opposed to galaxies with R_e ≈ 1 kpc that have Δv_Lyα ≈ 300 km s⁻¹. However, r_P = 0.48 and r_S = 0.58 for these data so at best the correlation is only tentative.

Finally, to judge if any systemic offsets of the resonant C IV emission was detectable with respect to the estimated systemic redshift, we looked at the velocity offsets of C IV for all C III emitters. Eleven of the thirteen z > 1.5 C III emitters also showing C IV show estimated Δv_C IV ≲ 250 km s⁻¹ within the error bars with a median of 92 km s⁻¹ (including the two outliers the mean is 82 km s⁻¹). We have avoided an attempt to disentangle the stellar and nebular contribution to the C IV emission (see Sect. 5), but these velocity offsets are in agreement with the stellar absorption based Δv_C IV estimates at lower redshift from Du et al. (2016). Hence, the estimated Δv_C IV with respect to systemic appears to be rather modest in the C III−C IV emitter sample. It is therefore also likely that the large scatter of the C IV measurements seen in the upper right panel of Fig. 11 can be contributed to mostly the Lyα velocity offsets when larger than ≈250 km s⁻¹. The two objects with Δv_C IV > 250 km s⁻¹ (720320277 and 720830605) show larger offsets due to prominent C IV P Cygni profiles as described in Sect. 5. Therefore, the estimated velocity offsets for these objects reflect the “residual” emission in the red part of the P Cygni profile. As mentioned, full modeling of the C IV line profiles is beyond the scope of this work.

9. Estimating electron density

The relative strengths of the C III and Si III emission line doublet components are highly sensitive to the electron number density of the emitting gas and the ISM pressure (e.g., Keenan et al. 1992; Osterbrock & Ferland 2006; Kewley et al. 2019a,b). Using the PyNeb software (Luridiana et al. 2013, 2015) we estimated the electron density, n_e, for the 52 and 13 objects with detected C III and Si III emission, respectively (see Fig. B.1). The doublet component flux ratios are determined directly from the FELIS template matches. In 0/13 and 7/52 of the Si III and C III doublet detections the fainter component is only marginally detected. However, for these few sources the flux ratio, and hence an estimate of the electron density, can still be inferred from the template matches, as the flux of each doublet component as mentioned varies independently. In these cases the flux ratio will naturally be pushed to the edges of the allowed ranges for the templates (see Table 2). The mean S/N(FELIS) of the doublet detections used to infer the electron density is 5.05 and 7.67 for the Si III and C III detections, respectively. The S/N values span the ranges 3.07−10.30 and 3.03−28.58 for the two samples. Example spectra and the corresponding FELIS template matches for a subsample of these sources are shown in Figs. 3, E.2, E.5, and E.6. Figure 16 shows the distribution of the n_e estimates from our MUSE samples assuming a fixed electron temperature of 10⁴ K. Fixing the electron temperature to 5000 K or 2 × 10⁴ K instead would shift the data points onto the lighter or darker shaded red curves shown in Fig. 16, respectively. The estimates sample the full range of electron densities where the C III and Si III ratios are sensitive, but the majority of the estimates (and upper limits) have n_e < 10⁵ cm⁻³. We checked for correlations between the inferred ISM electron densities and the estimated EW₀(C III) (and EW₀(Si III)) values, as larger EW₀(C III) potentially indicates younger systems, but found no clear dependencies.

Fig. 16. Electron densities (ne) estimated based on the C III (top) and Si III (bottom) doublet line flux ratios for a fixed electron temperature of Te = 104 K. The error bars (and limits) correspond to the 1σ uncertainties on the doublet flux ratios. For comparison, the theoretical curves (e.g., Osterbrock & Ferland 2006; Kewley et al. 2019a,b) for electron temperatures of 5000 K, 10 000 K and 20 000 K are shown as the solid red curves. The filled histograms count the number of objects where ne could be determined from the observed flux ratio. The open histograms add the subset of objects where the ne estimates are upper limits. The single open symbol to the right in the top panel indicates the C III flux ratio of object 605172634 (Fig. E.3) which results in an unconstrained electron density estimate for Te = 104 K.

Figure 17 shows a comparison of the electron density estimates for the objects with measurements from both C III and Si III (11 sources all from the UDF mosaic and UDF10). Within 3σ the estimated electron densities generally agree, even though several objects only have upper limits on n_e(C III) and/or n_e(Si III). These 11 sources have mean S/N(FELIS) of 5.34 and 14.52 spanning the ranges 3.07−10.30 and 5.35−28.58 for the Si III and C III detections, respectively.

Fig. 17. Comparison of electron density estimates based on the C III (x-axis) and Si III (y-axis) from Fig. 16 for the 11 objects (all from the UDF mosaic and UDF10) with measurements from both doublets. The shown electron densities were estimated for a fixed electron temperature of Te = 104 K. The errorbars (and limits) corresponds to the 1σ uncertainties on the doublet flux ratios.

The upper limits on n_e result from doublet flux ratios above the theoretically allowed values (solid red curves in Fig. 16) within the 1σ error bars on the doublet component flux ratios. The upper limits are then quoted as the electron density corresponding to the lowest value in the 1σ uncertainty range on the flux ratio as shown by the vertical error bars on the upper limits in Fig. 16. In some cases, if not just a result of large uncertainty on the measured flux ratios, these upper limits could indicate that the assumptions for the PyNeb calculations are incorrect. For instance, the outlier in the upper right corner of the top panel in Fig. 16 (right-hand side of Fig. 17) is object 605172634 (open symbol). This object (manually put at n_e = 7 × 10⁶ cm⁻³) has a well-constrained C III doublet component flux ratio of 1.7 ± 0.1 from the emission line template fit to the C III emission (shown in Fig. E.3). This leads to an un-constrained estimate of n_e as the 1σ error range is outside the allowed theoretical range for n_e given the assumption on the electron temperature. A higher electron temperature (darker red curve) allows for larger flux ratios. However, for a C III flux ratio of 1.7 ± 0.1 a very high electron temperature of roughly 10⁵ K is needed to get a constraint on the electron density (n_e < 5 × 10⁵ cm⁻³). In a similar fashion, some of the lower inferred UV emission line doublet flux ratios with large uncertainties quoted as upper limits in Fig. 16 can also be turned into estimates of an allowed electron density range by considering higher electron temperatures (and vice versa). A caveat of the classical density diagnostics that could also play an important role in explaining similar measurements is the existence of density (and temperature) inhomogeneities that challenges the simple interpretation of constant temperature and density throughout the emitting gas (Peimbert et al. 2017).

The electron densities estimated for the sample studied here are generally larger than the electron densities (obtained from [S II] λ6716 + [S II] λ6730) of the “green peas” and “Lyman break analogs” studied by Jiang et al. (2019), where they find n_e ≲ 10³ cm⁻³. However, as shown by the curves in Fig. 16 both C III and Si III saturate and are not sensitive to these low densities. At the same time the [S II] flux ratio saturates at n_e ≈ 10⁴ cm⁻³ (Kewley et al. 2019a) providing little overlap between the UV and optical tracers. The electron densities estimated via [S II] λ6716 + [S II] λ6730 and [O II] λ3726 + [O II] λ3729 for the sample of more common z ≈ 2.3 galaxies described by Sanders et al. (2016b) show n_e ≲ 3 × 10³ cm⁻³. Similar estimates were presented at z = 0.1−1 by Ly et al. (2016). Again, the optical probes used in these studies are not sensitive to n_e ≳ 10⁴ cm⁻³ (Kewley et al. 2019a). Hence, there appear to be a few objects with relatively high densities in our sample (similar to what is seen in the objects presented by James et al. 2014, 2018), whereas the large number of upper limits for the Si III and C III estimates are fully consistent with electron densities of special but also more generic galaxy samples at lower redshifts. However, as discussed by Maseda et al. (2017), considering the differences between electron density estimates obtained from lines with different ionization energies, one should keep in mind that these lines originate in physically different parts of the nebulae so different densities are to be expected even for the same parent galaxy. For example, the S II lines are generated in the outskirts, that is lowest density parts assuming a density stratification, of the ionized regions, while C III and Si III originate in the inner denser parts.

10. Estimating gas-phase abundances

The recent study by Byler et al. (2020) tested the diagnostics presented by Byler et al. (2018) using a sample of galaxies at redshift below 0.1 and at z = 2−3. They used predictions of UV emission fluxes from the Flexible Stellar Population Synthesis (FSPS) nebular emission model (Conroy et al. 2009; Conroy & Gunn 2010) to infer the connection between UV emission strengths and the gas-phase abundance parameterized as 12 + log₁₀(O/H). They compress the multidimensional space of predicted UV line fluxes from Si III, O III, He II, and C III into a set of fitting formulas that provide estimates of the gas-phase abundance. The first is based on a combination of Si III, O III, and C III dubbed Si3-O3C3:

(4)

where x corresponds to log₁₀(O III]₂/C III) and y is log₁₀([Si III]/C III). Byler et al. (2020) quote a typical statistical error for this relation of ±0.14 dex. The second fitting formula is based on a combination of O III, He II, and C III dubbed He2-O3C3:

(5)

where x is again log₁₀(O III]₂/C III) and y corresponds to log₁₀(He II/C III). The quoted typical statistical error for this fit is ±0.08 dex. We estimate the gas phase abundance following the Si3-O3C3 (He2-03C3) method for all seven (four) MUSE objects and 76 (82) objects from the collection of data from the literature described in Appendix C which have all relevant lines detected above a S/N of 3. Figure 18 shows the results applying the Byler et al. (2020) fitting formulas to these sources. The MUSE objects from our study are shown with large symbols; z ≈ 2−3 and literature sample as small dots.

Fig. 18. Gas-phase abundances estimated using the third order polynomial fitting functions for the Si3-O3C3 (top panel) and He2-03C3 (bottom panel) diagnostics from Byler et al. (2020) as a function of source redshift. The large symbols are MUSE sources from the UDF and small dots are from the literature collection described in Appendix C. The solid histograms show the subset of the objects from this work, whereas the dotted histograms also include the sources from the literature. All points are color coded according to EW0(C III). The gray bands mark regions of super-solar gas-phase abundances at 12 + log10(O/H) > 8.69.

First, we see that all estimates of the gas-phase abundances from our study and the literature sample show subsolar abundances. This indicates that sources with prominent UV emission lines are mostly systems with low gas-phase abundances, that is with low gas-phase metallicities. This agrees with the notion that the emitting gas of star forming galaxies capable of producing strong rest-frame UV emission generally has low (subsolar) metallicity with higher temperatures and stronger ionization fields of gas surrounding young massive stars (e.g., Stark et al. 2014; Maseda et al. 2017; Senchyna et al. 2019; Feltre et al. 2020). The fact that none of the UV emission line sources from the large collection of sources assembled from the literature show super-solar gas-phase abundances supports this. There is no indication that the z > 1 data would be particularly biased and only contain systems with low gas-phase abundances. In several cases the line emitters were detected in broader representative samples of galaxies. In the case of the MUSE sources studied here, the parent sample consists of line-emitters (mostly LAEs) from an unbiased untargeted emission line search in the full data cubes.

The z ≈ 0 observations from the literature do on average have slightly higher gas-phase abundances and EW₀(C III) than the measurements at higher redshift. However, as the EW₀(C III) color coding in Fig. 18 shows, these galaxies predominantly have strong C III emission which might cause a selection bias rather than illustrating an underlying intrinsic correlation between EW₀(C III) and 12 + log₁₀(O/H). This also results in an apparent slight decrease in gas-phase abundance with redshift, which is however fully consistent with a flat nonevolving gas-phase abundance as a function of redshift within the scatter and errors of the measurements. From gas phase abundance measurements based on rest-frame optical emission lines, it has been suggested that there is a deficiency of EW₀(C III) ≳ 10 Å objects at 12 + log₁₀(O/H) ≳ 8.5 (Maseda et al. 2017; Senchyna et al. 2017). Senchyna et al. (2021) argues that this deficiency might even be for objects at 12 + log₁₀(O/H) ≳ 8.0 (Z/Z_⊙ ≳ 0.2). In line with these findings the UV-based gas-phase abundances presented in Fig. 18, which are all below 8.5, do not present any clear evidence for a correlation between EW₀(C III) and 12 + log₁₀(O/H) as just explained. It is however worth noting that a comparison by Rigby et al. (2021) between optical and the Byler et al. (2020) UV gas phase abundance diagnostics for a single gravitationally lensed source find that the UV estimates are 0.5−0.8 dex lower than the optical estimates. A correction of this size to the measurements presented in Fig. 18 would shift the values to roughly 8.5, which is the scale where trends with EW₀(C III) starts to emerge when considering optical gas phase abundance diagnostics.

Figure 19 presents a direct comparison of the four MUSE sources and 57 objects from the literature where both estimates could be performed. All EW₀ and flux measurements for these subsamples of sources with close to the full suite of rest-frame UV emission lines detected are available for the MUSE objects in the catalog described in Table B.1 and for the literature sample in the catalog described in Table C.1. In Fig. 19 we see that the He2-03C3 gas-phase abundance estimator predicts lower values than the Si3-O3C3 estimator for Si3-O3C3 abundances of 12 + log₁₀(O/H) ≳ 8. This confirms the results from the similar comparison performed by Byler et al. (2020, their Fig. 11). Also the average gas-phase abundance predicted for the collected source samples for the two predictors (12 + log₁₀(O/H) ≈ 7.6 and ≈8.2 for He2-03C3 and Si3-O3C3, respectively) agrees with the estimates by Byler et al. (2020). We caution that it has been pointed out that the He II emission line is potentially a problematic tracer of the gas-phase abundances, as it likely includes both nebular emission and emission from stellar winds as described by Byler et al. (2018, 2020) among others. They however stress that He2-03C3 is a reliable metallicity tracer, particularly at low gas-phase abundances (12 + log₁₀(O/H) ≲ 8), where stellar contributions are minimal. Byler et al. (2020) furthermore find that the ability for He2-03C3 to reproduce the abundance estimates based on optical lines depends on the photoionization grid used for individual sources, which also indicates that indeed multiple radiation processes contribute to the He II emission.

Fig. 19. Comparison of the gas-phase abundance estimates from Fig. 18 for objects with all relevant emission features detected. The large symbols are MUSE sources from the UDF and small dots are objects from the literature collection described in Appendix C. The solid histograms show the subset of objects from this work and the dotted histograms include the sample from the literature. The gray region marks values of super-solar abundances (12 + log10(O/H) > 8.69). Points are color-coded according to object redshift.

11. Physical parameter inference from photoionization models

Photoionization models are not only useful for probing the gas-phase abundances as just described and exemplified. They also present predictions for a range of physical galaxy properties given theoretically predicted emission line fluxes. Thus, by constraining the emission line fluxes and ratios between individual line species observationally, model comparisons can provide constraints on the emitting galaxies and their environment’s likely physical properties. To infer the characteristics of the observed galaxy sample from MUSE studied here, we consider three suites of photoionization models, taken from the literature, that reproduce the nebular emission from different ionizing sources (massive single and binary stars and AGN) obtained combining the ionizing radiation field of these sources with the CLOUDY photoionization code (c13.03; Ferland et al. 2013). Specifically, the models we consider here are spectral models of the nebular emission from gas ionized by single young and massive stars by Gutkin et al. (2016), models that include the contribution from binary stars by Xiao et al. (2018), and the Feltre et al. (2016) models which describe the emission of the gas ionized by an AGN. In the following we summarize the main features of these models relevant for our work. For further details we refer to the individual papers describing the models.

The Gutkin et al. (2016) models, developed as part of the NEw frOntiers in Galaxy spectrAl modeLing (NEOGAL⁵) project, compute the nebular emission from the gas in H II regions using the latest version of the stellar population synthesis models of Bruzual & Charlot (2003) to be presented by Charlot & Bruzual (in prep.). These models incorporate updated stellar evolutionary tracks (Bressan et al. 2012), including new prescriptions for the evolution of the most massive stars in the Wolf-Rayet phase. The Gutkin et al. (2016) models provide the nebular emission from a whole galaxy, parametrized in terms of “galaxy-wide” parameters, by convolving the spectral evolution of single, ionization bound H II regions with a constant star formation history. Here we use the publicly available photoionization models, each computed at 90 stellar population ages between 0 and 10 Myr, assuming a constant star formation rate for 100 Myr, resulting in this being the effective age of the most recent episode of star formation in the models (Gutkin et al. 2016).

To explore the effect of binary interactions, we consider the Xiao et al. (2018) photoionization models, which provide the emission of H II regions ionized by single or binary stellar populations treated as a single instantaneous starburst with a given age. These models are based on the Binary Population and Spectral Synthesis (BPASS⁶) v2.1 models (Stanway et al. 2016; Eldridge et al. 2017; Stanway & Eldridge 2018, the latter referring to v2.2). We consider the available models for binary stellar populations of ages from 10⁶ to 10⁸ years. By including the BPASS-based models we attempt to accommodate the influence and likely important effect binary stellar populations have on the amount and strength of the ionizing photons being produced by star formation. In particular, for high-ionization lines like C IV and He II it has been argued that binary stellar populations are needed, though not always enough, to produce the observed emission from non-AGN galaxies (e.g., Steidel et al. 2016, 2018; Berg et al. 2019a; Nanayakkara et al. 2019; Senchyna et al. 2021).

To also account for the potential contribution of an AGN, we consider the emission from the narrow-line gas emitting regions in AGN computed by Feltre et al. (2016) (also part of NEOGAL), where the ionizing radiation field is represented by a broken power law, with the UV spectral index α for F_ν ∝ ν^α between five and 1000 eV.

Each of the models are parametrized in terms of physical quantities like ionization parameter (log₁₀U) at the edge of the Strömgren sphere (i.e., the dimensionless ratio between the number of ionizing photons and the total number of hydrogen atoms computed at the Strömgren radius), the hydrogen density of the gas cloud (log₁₀n_H), the carbon-to-oxygen abundance ratio (C/O; only sampled in the models by Gutkin et al. 2016), and the metallicity (Z), that is the mass fraction of all elements heavier than helium. Hence, this metallicity estimate, which therefore includes both gas-phase and interstellar (dust-phase) metals, differs from the gas-phase abundance metallicity (12 + log₁₀(O/H)) discussed earlier.

The NEOGAL⁷ (Gutkin et al. 2016; Feltre et al. 2016) models also include dust physics (e.g., van Hoof et al. 2004, for grain physics in CLOUDY) and a self-consistent treatment of metal abundances and dust depletion. The depletion of refractory metals onto dust grains is parametrized by means of the dust-to-metal mass ratio (ξ_d). We note that the treatment of dust and, therefore, of all the related effects, like photon absorption and scattering, radiation pressure, collisional cooling, photo-electric heating of the gas, and metal depletion, are not included in the Xiao et al. (2018) BPASS models (see, e.g., Shields & Kennicutt 1995; Dopita et al. 2002, for the impact of dust on the emergent nebular emission). A subgrid of models of the nebular emission from gas ionized by binary stars (using the BPASSv2.2) which includes dust physics and depletion is presented by Plat et al. (2019). As an overview of the various parameters used for the three suites of models considered here and to ease comparison between them, we summarize the sampled parameter spaces and their ranges in Table 5. This table summarizes the grid of models that predict the line intensities of galaxies and AGN based on stellar populations with both single and binary stars. As we searched the MUSE objects for N V, C IV, He II, O III, Si III, and C III and therefore have constraints on all of these lines (when they fall in the observed rest-frame wavelength range), in principle the predictions of the photoionization models span a fifteen-dimensional space of possible flux ratios consisting of the line ratios N V/C IV, N V/He II, N V/O III, N V/Si III, N V/C III, C IV/He II, C IV/O III, C IV/Si III, C IV/C III, He II/O III, He II/Si III, He II/C III, O III/Si III, O III/C III, and Si III/C III. Here we have combined the fluxes of the individual components of the UV emission line doublets, for instance, F(C III) = F([C III]) + F(C III]), to limit the number of dimensions spanned.

Being projections of the multidimensional model grid, 2D explorations and projections of line ratios do not provide a robust way of inferring the physical parameters of the objects considered and fail to convey the full information and constraints available for them. This can to some extent be remedied by “marginalizing” over individual parameters to lower the dimensionality of the model parameter space explored, for instance, by assuming a fixed carbon-to-oxygen ratio, hydrogen number density, or ionization parameter of the models to display. However, to condense the full information of the photoionization models we here present a new approach, where we explore the distributions of the model parameters of the full multidimensional photoionization model grid points (see Table 5) and their flux ratio predictions. We assume flat priors on the individual parameters for all galaxies, and let the data tell us what the most likely best-fit models are. The sample of best-fit models agreeing within, for instance, 3σ of the measured emission line fluxes (limits) provide a distribution of the model parameters sampled by this subset of models. We refer to these parameter distributions as PIM-PDFs (PhotoIonization Model Probability Density Functions) indicating that they resemble the “probability density functions” of the actual physical parameters of the individual objects given the measured constraints on the observed UV emission line flux ratios. Two examples of PIM-PDFs are shown in Fig. 20 (for reference and comparison Appendix D shows the “full” PIM-PDFs for the NEOGAL and BPASS-based photoionization models in Fig. D.1). The PIM-PDF approach is different from the inferences performed by codes like the BayEsian Analysis of GaLaxy sEds (BEAGLE, Chevallard & Charlot 2016) and Prospector (Leja et al. 2017; Johnson et al. 2021). These codes sample the true PDF through Bayesian inference and Markov chain Monte Carlo for each object analyzed. Hence, these codes are optimizing the parameter grid as opposed to the PIM-PDFs which by construction just select among the available precomputed models and base the resulting distribution on these estimates. This approach is what makes PIM-PDFs more efficient when evaluating large samples of objects, as each object is not independently forced through a time-consuming optimization process.

Fig. 20. Examples of PhotoIonization Model Probability Density Functions (PIM-PDFs) for the objects 102014087 (top) and 601281436 (bottom). The histograms in the different panels show the distribution of AGN (purple; Feltre et al. 2016), star forming (blue; Gutkin et al. 2016), BPASS binary (green) and single star (yellow; Xiao et al. 2018) photoionization model parameters satisfying the observational emission line flux ratio constraints listed in the left margin. The spacing of the parameter histograms’ bars indicate the discreteness of the sampling of each parameter space sampled by the models as listed in Table 5. The spectra of the two objects are shown in Fig. E.3. The limits on the observed flux ratios are all 3σ. The vertical dotted lines show the median values and the 68% confidence intervals of the distributions. The vertical dashed lines mark the mean values. For comparison the full unconstrained set of PIM-PDFs are shown in Fig. D.1.

The first object (102014087; top panel) is a case where the observational constraints (listed in the left margin) are more easily reproduced by models of ionizing photons from star formation as only 0.08% of the available AGN models can reproduce the observed emission line flux ratios. The set of models including contribution from binary stars (green PIM-PDFs) have a higher fraction of solutions that are capable of reproducing the observations. The resulting PIM-PDFs predict that the object 102014087 likely has a (relative) metallicity of roughly a few times 10⁻² (Z_⊙ = 0.01524; Gutkin et al. 2016; Feltre et al. 2016) and an ionization parameter that is ≳10⁻³. The stellar populations of the BPASS models prefer ages above ≈10⁷ years. The hydrogen number density and the additional parameters sampled by the NEOGAL model grids are essentially unconstrained for the models that are able to reproduce the observations of 102014087, even though the preferred C/O ratio appears to be small rather than large.

The observational constraints for the second object shown in Fig. 20 (601281436; bottom panel) have a larger set of models with an AGN as ionizing source capable of reproducing the observational constraints, as essentially all the NEOGAL AGN models (purple) can reproduce the observations (and therefore poorly constrain the physical parameters), whereas the NEOGAL star formation model grids (blue) are struggling as only 1.38% of them can reproduce the measured flux ratios. The few models that can have low C/O and a high mass-function cutoff. A subsample of roughly 10% of the BPASS-based models with preferentially small n_H and age, a relative metallicity of 10⁻⁵ and log₁₀(U) ≲ −2.5 are also capable of reproducing the observational constraints. The spectra of both 102014087 and 601281436 are shown in Fig. E.3.

We produced the PIM-PDFs for all studied objects with flux ratio constraints, that is the 103 objects with at least one UV emission line detected. Among these objects 62 (for 23 of these sources the line has S/N(FELIS) > 5) only have a single line detected, resulting in an inference based on flux ratio limits only. By recording the mean values and the standard deviation of the PIM-PDFs for all objects in our sample, we can estimate the distribution of the most likely physical (model) parameters for the objects detected in our study. Figure 21 presents these statistics for the ionization parameter (log₁₀(U), left panel) and metallicity (Z, right panel) for the full sample. Points at the top are generally unconstrained given the large relative width of the PIM-PDFs (σ/|mean|), whereas points towards the lower parts of each panel provide reliable estimates of the sampled parameter. For the distributions of the ionization parameter in the left panel the preferred values appear to scatter around log₁₀(U)≈ − 2.5. In the right panel we see that the full sample of PIM-PDFs predicts systems (of star-forming models) that span a range in metallicity from roughly 5 × 10⁻³ to 3 × 10⁻² (corresponding to 0.32 ≲ Z/Z_⊙ ≲ 1.97) with a mean value of approximately solar at Z_⊙ = 0.01524 indicated by the vertical gray dashed line. As mentioned earlier, Tang et al. (2021a) argue that EW(C III) > 10 Å mostly requires log₁₀(U)≲ − 2.5 and metallicity of 0.2 Z_⊙ or lower. If we add the results from the PIM-PDFs based on the emission line fluxes from the collection of UV line emitters from the literature, the distribution of the best-fit log₁₀(U) stays roughly the same, whereas the Z distribution shifts to slightly lower values with a mean of subsolar metallicity at Z ≈ 10⁻² (Z/Z_⊙ ≈ 0.7). Hence, the PIM-PDFs predict the larger sample of MUSE galaxies studied here to have higher metallicities compared to metallicity estimates of samples at similar redshifts in the literature (e.g., Maiolino et al. 2008; Sanders et al. 2020).

Fig. 21. Mean value (x-axis) and normalized standard deviation (y-axis) of the ionization parameter (log10(U)) and the metallicity (Z) PIM-PDFs for the MUSE-Wide, the UDF mosaic, and the UDF10 objects studied in this work. The blue points represent the distribution of star forming models (Gutkin et al. 2016) matching the observations, the purple points indicate AGN models (Feltre et al. 2016), the green points show models based on BPASS binaries (Xiao et al. 2018), and the yellow points represent models of BPASS single stars (Xiao et al. 2018). The solid histograms show the distribution of the individual parameters. The vertical dashed line in the right panel marks the solar metallicity Z⊙ = 0.01524.

We performed similar comparisons for the remaining parameters sampled by the photoionization models listed in Table 5. The MUSE sample as well as the sample from the literature (Appendix C) prefer values of m_up = 200 M_⊙, ξ_d ≈ 0.3 and α ≈ −1.6 from the best-fit NEOGAL models. The BPASS-based model’s best-fit age distributions peak at the central value of 10⁷ years for both the MUSE and literature sample. The C/O abundance ratio from the Gutkin et al. (2016) models spans the full range of sampled values for the MUSE objects with a peak around the solar value (C/O)_⊙ = 0.44. The sample of emitters from the literature more strongly favors C/O ratios around solar. The distributions of the neutral hydrogen number density for the combined sample PIM-PDFs of the BPASS and NEOGAL star formation models peak at values between 10² cm⁻³ and 10³ cm⁻³. The solutions from the BPASS-based models generally predict densities lower than the NEOGAL models, which reflects the different sampling ranges of the two different suites of models as the means of the sampled ranges are 10^1.5 cm⁻³ and 10³ cm⁻³ for the BPASS and NEOGAL models, respectively.

Even though the inferences on the physical parameters of the individual objects and the galaxy samples using the PIM-PDFs presented here estimates the characteristics of the UV line emitters, it should be stressed that these predictions should be considered only indicative of the true intrinsic values. First, there are uncertainties and assumptions in producing the photoionization models resulting from the sampled parameter grids. Second, the recent study by Florian et al. (2021) presents spatially resolved estimates of physical parameters based on comparisons between rest-frame optical lines and photoionization grids of two lensed (and hence highly resolved) sources from HST grism spectroscopy. They show and discuss the challenges in inferring metallicity and ionization parameters based on such grids, both in terms of uncertainties in the models, but also in terms of intra-object differences between, for instance, individual star-forming regions and integrated estimates of line ratios. This of course also applies to the inferences on the photoionization model parameters obtained for the objects studied here, where all quantities are integrated over the whole galaxy given the nonlensed nature of the objects.

Another limitation of this approach is the nonuniform and in some cases coarse sampling of the physical parameters considered. This is done to limit the number of models to generate when calculating the emission line flux predictions when producing the photoionization model outputs. However, this means that well-constrained line fluxes (with small uncertainties) translate into small regions of parameter space being allowed by observations. In those cases, it might very well be that no models are able to reproduce the observations to within 3σ (as used here) of the measured line flux ratios. This does not mean that the models fail, but rather that they are not sampled finely enough to reflect the precision of the measurement. This issue of course scales with the number of well-constrained lines. A way to solve this is to either increase the considered multidimensional region of parameter space used to generate the PIM-PDFs, produce new more finely sampled models, or consider only the single best-fit model (in terms of χ²). Alternatively, dedicated machinery to infer physical parameters of individual objects, like for instance BEAGLE, can be used. A detailed comparison between the PIM-PDF approach presented here and other methods, is beyond the scope of the current work, but will likely prove useful and insightful when exploring further the physical properties of individual galaxies from the MUSE and literature samples in the future.

Lastly, more fundamental challenges with the considered photoionization models also provide empty PIM-PDFs when attempting to match observational constraints. An example of one such challenge or limitation to the models in general, and therefore also to the PIM-PDF approach, is He II emitted in star-forming galaxies. As we noted in Sect. 10, He II emission can contain contributions from both nebular emission and stellar winds (e.g., Shirazi & Brinchmann 2012; Nanayakkara et al. 2019; Byler et al. 2018, 2020), and producing enough high-energy photons to produce the He II emission has been proven challenging. Binaries (for example, as implemented in the BPASS models) have in some cases relieved part of the tension (Steidel et al. 2016), but they are often not enough. It has therefore been proposed and debated that the inclusion of X-ray binaries (Schaerer et al. 2019; Senchyna et al. 2020) or energetic shocks (Allen et al. 2008; Senchyna et al. 2017; Plat et al. 2019) could provide means for reproducing the observed He II fluxes and other lines produced by high-energy ionizing photons like N V and C IV within photoionization model frameworks.

Despite these limitations and challenges with inferring physical properties from photoionization models based on a set of observation constraints, the PIM-PDFs presented here provide a simple and fast way to obtain photoionization model predictions for large samples of objects including all available information and avoiding marginalization over individual parameters or projections onto lower-dimensional parameter spaces.

12. Conclusions

We have presented a comprehensive search for rest-frame UV emission lines in a parent sample of 2052 emission line galaxies detected with the MUSE integral field spectrograph in the COSMOS, GOODS-South, and HUDF as part of the MUSE-Wide and MUSE-Deep GTO surveys (Sect. 2). The studied objects were selected to have redshifts above 1.5 based on the emission feature(s) identified during the untargeted emission line search in the MUSE data cubes. This resulted in 1997 LAEs at z > 2.9 and 55 objects in the so-called MUSE redshift desert at 1.5 < z < 2.9 where neither [O II] λ3726 + [O II] λ3729 nor Lyα are available for source redshift identification (Sect. 3). For each of the objects in the parent sample, we extracted 1D spectra from the 3D MUSE data cubes optimized in both flux and S/N using the software TDOSE (Sect. 4). These spectra were searched for rest-frame UV emission lines red-wards of Lyα using Gaussian emission line template matching using the tool FELIS (presented as part of this paper in Sect. 5 and Appendix A). We visually vetted all potential detections of the UV emission lines N V, C IV, He II, O III, Si III, and C III and found 54 line emitters with 3σ detections among the 1997 LAEs. Including the sample in the redshift desert this number increases to 103 line emitters. Table 3 summarizes the complete sample of rest-frame UV detections from our search which are made publicly available with this paper (Appendix B). To further improve the statistics of the assessments and for comparison purposes, we complemented the main sample with an extensive collection of UV emission line flux measurements from the literature which are also made available with this work (Appendix C).

Based on the UV emission lines and complementary measurements of LAE characteristics from Kerutt et al. (2021), we explored the range of physical parameters of the galaxies in our main sample. Our main conclusions can be summarized as follows:

• The fraction of objects with detected rest-frame UV emission lines grows with increasing depth of the data. This implies a relative change in the shapes (slopes) of the luminosity function of the parent sample and the luminosity function of the UV line emitter subsamples (Sect. 5).

• The strengths of He II, O III, and Si III correlate with the flux of the C III emission as parametrized by the correlations presented in Table 4.

• Based on the FELIS Gaussian template matches we determine the relative strength of the UV emission line doublet component flux ratios. Here we find that F(O III]₁)/F(O III]₂), F([Si III])/F(Si III]), and F([C III])/F(C III]) are all in agreement with expected theoretical values for the vast majority of objects (Sect. 5 and Table 4).

• We determine EW₀ for all identified UV emission lines, and similar to the flux measurements we find significant correlations between them for the He II, O III, Si III, and C III emission lines (Sect. 6).

• We find correlations between EW₀(He II), EW₀(Si III), and EW₀(C III) and EW₀(Lyα). The majority of the LAEs with C III detected have EW₀(C III) ≈ 0.22 ± 0.18 × EW₀(Lyα) for 10 ≲ EW₀(Lyα)/Å ≲ 100 which is in agreement with previous studies (Sect. 7).

• Considering the subsample of LAEs only, we explored correlations between multiple LAE characteristics including EW₀(Lyα), M_UV, spectral slope, and LAE effective radius and found no prominent relationships. We do however see potential correlations between Lyα luminosity and the UV emission lines indicating that objects with lower Lyα luminosity have larger UV emission line EW₀ estimates (Sect. 7).

• The detection of multiple rest-frame emission lines enabled us to assess the velocity offset Δv of resonant lines like Lyα and C IV with respect to the systemic redshift probed by, for instance, C III (Sect. 8). In agreement with previous measurements we find Δv_Lyα of 250−500 km s⁻¹. Again checking for correlations with LAE characteristics we find broad agreement between theoretical and empirical relations between Δv_Lyα, M_UV, and EW₀(Lyα), despite a large scatter in the parameters for the objects studied here. In addition we confirm the empirical relation between the Lyα line width, Lyα peak separation, and Δv_Lyα even though a large scatter in the individual correlations is present. For the resonant emission of C IV we find that Δv_C IV ≲ 250 km s⁻¹ with a few outliers caused by the limitations of trying to model C IV P Cygni profiles from combined absorption and emission by a pair of Gaussians.

• For objects with detected C III and Si III the estimated electron density from the doublet flux ratios is generally n_e < 10⁵ cm⁻³ (Sect. 9).

• Using the fitting formulas by Byler et al. (2020), we show that the sample of objects with simultaneous detections of He II, O III, Si III, and C III have subsolar gas-phase abundances of 12 + log₁₀(O/H) ≈ 8. We do not find any trends with redshift or EW₀(C III) for the estimated gas-phase abundances. In agreement with Byler et al. (2020) we find that the tracer including He II, O III, and C III in most cases results in lower metallicity estimates than the tracer based on O III, Si III, and C III (Sect. 10).

• Finally, we present a new approach to condense information from physical parameter grids of photoionization models taking the full amount of information into account without marginalizing over individual parameters. We refer to the resulting distribution of model parameters that are able to reproduce the observational constraints from the UV emission lines as PIM-PDFs (PhotoIonization Model Probability Density Functions; Sect. 11). The PIM-PDFs provide the distribution of the best-fit (integrated) model parameters of individual sources given the observational constraints on the UV emission line fluxes and flux ratios obtained from the FELIS template matches. We show that the general sample of emitters recovered from the MUSE data have an average ionization parameter log₁₀(U)≈ − 2.5 and a mean metallicity of order solar, though individual objects span the range 5 × 10⁻³ ≲ Z ≲ 3 × 10⁻² corresponding to roughly 0.32 ≲ Z/Z_⊙ ≲ 1.97 which is at the high end when comparing to previous estimates from the literature.

In summary, with the large sample of emission line sources recovered from the MUSE GTO surveys, we have demonstrated the wealth of information and physical properties that rest-frame UV emission features red-wards of Lyα probe. Apart from gaining insight into the individual galaxies and samples themselves, the rest-frame UV emission lines also provide promising probes of high-redshift galaxies at epochs where the Lyα emission is strongly affected and absorbed by the increasing neutral CGM and IGM. Especially in light of upcoming near-infrared missions like the James Webb Space Telescope (JWST) and the Nancy Grace Roman Space Telescope (formerly known as WFIRST), further exploration of the characteristics of rest-frame UV emission line sources can serve as benchmarks, links to lower redshift, and means of comparison and redshift evolution assessment.

Acknowledgments

We would like to thank Charlotte Mason for useful discussions and for providing the data for the curves shown in Fig. 13 and Dawn Erb for providing the observational data for the comparison sample studied by Steidel et al. (2014), also shown in Fig. 13. This work has been supported by the BMBF grant 05A14BAC and we acknowledge support by the Competitive Fund of the Leibniz Association through grant SAW-2015-AIP-2. AF acknowledges the support from grant PRIN MIUR2017-20173ML3WW_001. JS acknowledges the support from Vici grant 639.043.409 from the Dutch Research Council (NWO). GM received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No MARACAS – DLV-896778. This paper is based on observations collected at the European Organisation for Astronomical Research in the Southern Hemisphere under ESO programmes 094.A-0289(B), 095.A-0010(A), 096.A-0045(A), 096.A-0045(B), 094.A-0205, 095.A-0240, 096.A-0090, 097.A-0160, and 098.A-0017. This paper also makes use of observations made with the NASA/ESA Hubble Space Telescope obtained at STScI. This research made use of the following programs and open-source packages for Python and we are thankful to their developers: DS9 (Joye & Mandel 2003), Astropy (Astropy Collaboration 2013, 2018), APLpy (Robitaille & Bressert 2012), iPython (Pérez & Granger 2007), numpy (van der Walt et al. 2011), matplotlib (Hunter 2007), and SciPy (Jones et al. 2001).

References

Appendix A: Finding Emission Lines In Spectra (FELIS)

FELIS (Finding Emission Lines In Spectra) is a publicly available⁸ Python tool enabling the search for emission lines and other spectral features in 1D spectra. FELIS was developed specifically for our particular science project but offers a general tool suitable for general detection of (weak) features in 1D spectra. As described below, FELIS provides tools to build spectral templates and mock spectra (Appendix A.1) which can be matched to observed spectra providing S/N estimates of potential emission features (Appendix A.2). In Appendix A.3 we describe the results from testing FELIS on a set of idealized MUSE-Wide mock spectra with UV emission lines similar to the spectra analyzed in the study described in this paper.

The search for emission features with FELIS is performed via standard template matching by minimizing the χ² between the input data (spectrum) and the model (idealized template of spectral feature). FELIS provides an estimate of the significance of the model template match by providing the S/N for the minimized χ². Hence, FELIS matches the defined templates to the observed spectra by minimizing the χ² expression

(A.1)

Here the index i runs over the individual pixels of the spectrum, 𝒟 represents the data, that is the flux measurements of the spectrum, σ² is the variance on the pixel fluxes, 𝒯 indicates the template (normalized to an integrated flux of 1) to match to the spectrum, and α is the flux scaling to apply to the template to obtain the best match. Using that

(A.2)

(A.3)

for the minimum χ² value implies

(A.4)

with an estimated variance on α given by

(A.5)

Hence, the best match of the normalized template 𝒯 is obtained by scaling it by the flux α. The S/N of this template match can be estimated as

(A.6)

Cross-correlating each template with the spectrum, that is shifting the template over the spectrum pixel-by-pixel while minimizing χ² provides a χ² and S/N curve for each template and spectrum pair. Examples of these curves are shown in the two lower panels of Figure A.1). The maximum of the S/N curve indicates the optimal alignment of the template with the spectrum. Comparing the (S/N)_𝒯, max with the maximum S/N values for a range of templates matched to the same spectrum provides the overall best match to the measured fluxes in the spectrum, (S/N)_max, given the set of templates explored.

Fig. A.1. Overview of the FELIS template fit to a simulated spectrum of a C III emitter with a doublet flux ratio of F([C III])/F(C III]) = 1.4 corresponding to an electron density < 104 cm−3 (Osterbrock & Ferland 2006). The top panel shows the simulated spectrum (black) with the best-fit scaled template plotted on top (green). The second panel shows the S/N spectrum, that is flux/. The third panel from the top shows the template cross-correlation S/N calculated with Equation A.6. The maximum value is indicated by the vertical dashed red line in all four panels. The bottom panel shows the χ2 estimate of the template cross-correlation.

A.1. Building spectral templates with FELIS

The templates used for the FELIS cross-correlation can be generated with the function

felis_build_template.build_template()

The input to this function is a dictionary defining the template components to combine to the final template and the wavelength range and resolution of the template to generate. Various template components can be added to the dictionary including:

• A skewed Gaussian emission line profile of the form

  (A.7)

  (A.8)

  (A.9)

  where PDF and CDF refers to the probability and cumulative distribution functions, respectively. The template dictionary parameters to provide are

  ['GAUSS', mean, sigma, skew, scaling, 'info']
  

  where he parameters mean and sigma are represented by μ and σ. The β provides the skew of the Gaussian distribution. If β = 0 then the error function (erf) is 0 and G(λ) is a regular Gaussian PDF.

• A delta function (single pixel flux) added at the wavelength nearest the position in the generated template. The template dictionary parameters to provide are

  ['DELTA', position, total flux, 'info']
  

• A Gaussian line spread function (LSF) which the generated template will be convolved with before it is returned. The template dictionary parameters to provide for this are

  ['LSF', sigma, 'info']
  

  For the study presented in this paper, no LSF was added to the templates. The templates are “raw” Gaussians representing the observed signal focusing on recovering line fluxes.

• A linear continuum with the template dictionary parameters

  ['CONT', level, slope, lam0, 'info']
  

  describing the continuum as

  (A.10)

  where lam0 provides a reference wavelength λ₀ for the slope. Alternatively, a predefined continuum can be aded as a “flux feature”.

• A predefined spectral flux feature defined by a wavelength and flux vector as provided in the template dictionary parameter list:

  ['FEATURE',wavelength,flux,'info']
  

  Examples of such features could be spectral breaks, nongaussian emission line profiles, or a nonlinear continuum.

The “info” provided to each of the template components contains a string with information that will be written to the header of the fits file containing the template returned by FELIS. As an example generating a template containing a C III emission line doublet, a delta-function spike at 1900 Å, and a blue continuum can be done with the commands:

import felis_build_template as fbt
# --- INPUT ---
# fill template component dictionary
tcdic = {}
tcdic['D1900'] = ['DELTA', 1900.0, 10.0, 'Delta␣function␣
     ↪ at 1900A']
tcdic['CIII1'] = ['GAUSS', 1907.0, 0.5, 0.0, 10.0, 'CIII
     ↪ ]1907A']
tcdic['CIII2'] = ['GAUSS', 1909.0, 0.5, 0.5, 5.0, 'CIII
     ↪ ]1909A']
tcdic['CONT'] = ['CONT', 1.0, -0.03, 1908.0, 'Continuum␣
     ↪ with␣flux␣1.0␣at␣1908␣+␣slope␣-0.03']
# --- COMMAND ---
template_range = [1870,1980,0.1]
fbt.build_template(template_range,tcdic,tempfile='./
     ↪ outputname.fits')

This template will be noise-free. To enable the generation of mock spectra a noise description can be provided to

felis_build_template.build_template()

This adds noise to the template spectrum according to one of the following prescriptions:

• Noise drawn for each template pixel from a Poisson distribution around a mean value can be done with

  noise=['POISSON',mean]
  

• Gaussian random noise of each template pixel from a Gaussian error spectrum defined by the mean and standard deviation sigma in each pixel can be obtained with

  noise=['GAUSS',mean,sigma]
  

• Assigning a Gaussian random noise around a constant noise level with the value value to each template pixel is done with

  noise=['CONSTANT',value]
  

• A Gaussian random noise for each template pixel drawn from a provided error spectrum defined by the wavelength and flux vectors wave and flux after interpolating it to the template wavelength range can be obtained with

  noise=['SPECTRUM',wave,flux]
  

  This noise prescription can be used to include information about sky-residuals, poor pixels, etc. in the noise prescription.

In the last three noise prescriptions an error spectrum is defined (or explicitly provided) from which Gaussian noise is drawn at each pixel. In practice, the pixel values in the error spectrum (which are stored as the template’s flux error) are assumed to represent the standard deviation of the noise to be applied to each pixel. The noise-free template flux of each pixel is modified by adding a flux value drawn from a Gaussian distribution with mean 0 and standard deviation corresponding to the error spectrum value for that pixel. Hence, the resulting S/N of a featureless noisy template is ensured to be a Gaussian centered on 0 with a standard deviation of 1.

A.2. Searching spectra with FELIS templates

Having build a set of spectral templates that represents the parameter space of the spectral features expected to occur in the observed spectra, the spectra can be searched with FELIS by executing the following command:

import felis
# --- INPUT ---
spectra = [spectra,to,match]
redshifts = [1.5,3.5,5.5]
templates = [templates,to,match]
picklefile = './NameOfOutput.pkl'
plotdir = './path/to/directory/with/plots/'
wavewindow = [50]*len(spectra) # observed frame
windowcenter = [5007]*len(spectra) # rest-frame
# --- COMMAND ---
ccdic = felis.match_templates2specs(templates, spectra,
      ↪ redshifts, picklefile, plotdir=plotdir, wavewindow=
      ↪ wavewindow, wavecen_restframe=windowcenter)

This will produce a dictionary (ccdic) containing the results from the cross-correlation χ² minimizations performed by FELIS. The dictionary is saved to a binary Python pickle file if an output filename is provided. This pickle file can be loaded with

felis.load_picklefile()

To select spectra from the pickle file output that have template matches of a certain quality and characteristics the function

felis.selection_from_picklefile(picklefile)

is provided. For instance,

speclist = felis.selection_from_picklefile(picklefile,
     ↪ S2Nmaxrange=[3,5])

will return all spectra that have a template match with a maximum S/N between three and five. When provided a FELIS output dictionary and a dictionary key (the name of a spectrum) the function

felis.getresult4maxS2N()

will return information (including best-fit redshift, maximum S/N value, and template name) of the template matching the observed spectrum best. The content of the output can be plotted with

felis.plot_picklefilecontent()

Above, the redshifts are provided as input to perform the template matching by FELIS in rest-frame. Hence, the templates are assumed to be in rest-frame, and the observed spectra will be moved to rest-frame by scaling wavelength and fluxes such that:

(A.11)

(A.12)

(A.13)

This ensures that integrated fluxes and S/N are conserved between the observed-frame and rest-frame spectra.

To avoid cross-matching a template over the whole spectral range, if for instance the redshift of the object is roughly known and the search is for a well-defined emission line, the wavewindow and wavecen_restframe can be used to define the window over which the template match (cross-correlation) is performed. The window is defined as

(A.14)

The redshifts can be both photometric or spectroscopic. Supplying more uncertain photometric redshifts of course requires searching for features over a larger wavelength range. To look for potential velocity shifts of emission lines with respect to each other or systemic redshift, a spectroscopic redshift is clearly preferred.

A.3. Testing FELIS on MUSE-Wide mock spectra

To test the reliability of FELIS to recover line fluxes, line widths and doublet flux ratios at different emission line S/N we generated a series of mock spectra for regions around the main rest-frame UV emission lines studied in this work (see Table 2). The mock spectra were generated using the template creation capabilities of FELIS described in Appendix A.1. Each line was modeled as a Gaussian with a width σ_Gauss and a total emission line flux in units of 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹. For doublets the emission line flux ratios were set as an extra constraint. For the testing we also generate mock Lyα spectra where we represented the Lyα emission by a single skewed Gaussian ignoring any blue bumps. All spectra were generated with a wavelength resolution of 1.25 Å and generated in a wavelength range of ±50 Å around the central wavelength of the emission line (for doublets we used the central wavelength of the doublet). Table A.1 summarizes the template parameters for the mock spectra.

We generated a set of noise-free templates as control sample, and a sample of templates with idealized noise added. Each of the noisy mock spectra was added noise based on the median noise spectrum from MUSE-Wide (shown by the green and blue curves in Figure A.2) scaled by a factor 5.5. The scaling was introduced to resemble noise for a spectrum extracted using a aperture, which corresponds to propagating the noise from 30 MUSE spaxels as for N = 30. The flux value in each pixel of the mock spectra was obtained by adding the noise-free template a flux term obtained from drawing a random value from a gaussian centered around 0 with a width corresponding to the value from the MUSE-Wide noise spectrum at each individual wavelength.

Fig. A.2. Median noise spectrum of two random MUSE-Wide pointings (green) and 70 GOODS-South and COSMOS MUSE-Wide pointings (blue – almost identical to the green curve). The latter spectrum was used to draw noisy flux values for the MUSE-Wide mock spectra that FELIS was tested on. For comparison, the median noise spectra of the deeper MUSE UDF mosaic and the UDF10 pointing are shown in red and gray. The median noise spectra shown in this panel are displayed by the black filled regions at the bottom of each panel in the spectral overviews in Figure 3 and Figures E.1–E.6.

To recover the emission lines of the mock spectra we performed a search with FELIS in a ±10 Å window in the rest-frame around the central wavelength of the mock spectra of the targeted emission features. The mock spectra were positioned at arbitrary observed redshifts to include rest-frame conversions of line widths and interpolation of templates to the mock spectra rest-frame reference grids. We first matched the templates used in the main text described in Table 2 on the noise-free mock spectra. In this case the FELIS template matches recovers the input flux values (total flux of emission doublets) to well within 10% as long as the spectral feature has an integrated S/N above three and the intrinsic rest-frame line width of the feature in the mock spectrum is within the sampled range of the templates listed in Table 2. Also the recovered line widths, σ_Gauss, are accurate to within 10%. Note however, that the template resolution of σ_Gauss is 0.1, so for mock spectra with intrinsic rest-frame line widths below ≈0.5 Å the relative inaccuracy grows to roughly 20% for the noise free mock spectra. The intrinsic doublet flux ratios of the emission line doublets are recovered correctly to within a few percent (arising from the interpolation of the templates to the “observed” mock spectra’s rest-frame grid). In the recovery of the intrinsic parameters of the noise-free mock spectra, there is a correlation between the precision of the recovered total line flux and the intrinsic line widths such that for intrinsic rest-frame line widths comparable to the template resolution (here 0.1) the recovered flux is less precise.

Having fitted the mock spectra without noise we ran FELIS with the same set of templates (Table 2) on the same set of mock spectra (Table A.1) with noise added as described above. As expected the ability of the FELIS template match to recover the intrinsic flux values and emission line widths decreases when introducing noise in the mock spectra. Figure A.1 shows an example of one of the results from a template match to a noisy mock spectrum performed with FELIS. The mock spectrum (black curve) shows a C III-emitter with a doublet flux ratio of [C III]/C III] = 1.4 and a total intrinsic (before adding noise) line flux of 1028.57 × 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹. The noise in each pixel was drawn from a median effective noise spectrum of the MUSE-Wide fields at an “observed” redshift of 2.7 (i.e., at λ_obs ≈ 7055 Å). The intrinsic “observed” frame Gaussian line width was set to 2.0 Å which corresponds to a rest-frame line width of 0.54 Å. For the noisy mock spectrum the integrated line flux (integrated over [λ₁ − 3σ, λ₂ + 3σ] with σ = 2 Å) is 989.60 × 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹. The scaled best-fit template from the set of templates in Table 2 is shown in green. The best-fit template, returning the maximum cross-correlation (S/N)_𝒯, max ≈ 12 as shown in the third panel of Figure A.1, has a doublet flux ratio of 1.2, and a rest-frame line width of 0.5 Å. The total flux estimated by the template scaling performed by FELIS is (966.93 ± 80.82)×10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹. In the case of matching the templates from Table 2 to the noise-free version of the mock spectrum shown in Figure A.1 the estimates of rest-frame line width, flux ratio, and total flux were 0.5 Å, 1.4, and 993.47 × 10⁻²⁰ erg s⁻¹ cm⁻² Å⁻¹, respectively. The roughly 3% loss in total flux when comparing the noise-free template match to the intrinsic noise-free mock spectrum comes from the mismatch in rest-frame emission line widths, as the line widths of the templates take steps of 0.1 Å. Hence, only templates with 0.5 Å and 0.6 Å exist, and the best-fit template thus has a line width of only 0.5 Å compared to the actual rest-frame line width of the mock spectrum of 2.0 Å/(1+2.7) ≈ 0.54 Å. This is illustrated in the panels of Figure A.3 showing both the noise-free (top) and noisy (bottom) mock spectra with the best-fit scaled template shown on top in red.

In general, for the template match to the noisy mock spectra we see that the total flux is recovered with an accuracy better than 50% for the spectral doublets and to within 20% for single lines for a FELIS match significance S/N > 3. A FELIS S/N > 3 can roughly be translated to a total integrated flux of the spectral feature of ≈300−500 × 10⁻²⁰ erg s⁻¹ cm⁻² with the added noise from the median MUSE-Wide spectrum. This flux level is of course lowered for deeper data. The flux ratios are generally robust to within a factor of two of the intrinsic flux ratio, but are in most cases good at the ≈50% level, especially if a slightly higher S/N than three is required. The precision of the recovery of the flux ratio appears to be controlled by the intrinsic line widths and their interpolation to the rest-frame grid of the “observed” mock spectra, as the doublet separations are at least 4.8 times that of the rest-frame pixel resolution given the 1.25 Å spectral resolution of MUSE and the emission line doublets considered here. The Gaussian line width is recovered within 50% of the intrinsic value of the mock spectra. Again this assumes a FELIS S/N > 3 of the template match. Lastly, the redshifts, that is the location in wavelength of the best-fit template, are recovered within ±100 km s⁻¹. As no velocity offsets were included in the mock spectra, this shows that the expected precision of the recovered velocity offset from FELIS given the spectral sampling of MUSE (1.25 Å) is at the 100 km s⁻¹ level for the doublets considered here assuming a FELIS S/N > 3.

Fig. A.3. Noise-free (top) and noisy (bottom) template matches to the mock spectrum shown in Figure A.1. The top panel of Figure A.1 shows the same as the bottom panel of this figure. The red curve shows the FELIS template match to the part of the spectrum in black (within the vertical dotted lines). The green curve shows the template match before applying the flux scaling. The part of the mock spectrum outside the considered wavelength range is marked in gray.

Appendix B: Rest-frame UV emission line catalog and measurement correlations

Together with the paper we provide a catalog of the full list of MUSE emission line sources at 1.5 < z < 6.4 searched for UV emission lines in this paper. The content of this catalog is described in Table B.1.

As presented throughout the paper, several correlations between the UV emission line fluxes and EW₀ from this catalog (combined with the literature comparison sample described in Appendix C) were found. Among those are the correlations between the flux of the individual components of the emission line doublets shown in Figure B.1 and discussed in Section 5. The correlations between the UV emission line EW₀ estimates described in Section 6 are shown in Figure B.2. At the end of Section 7 we also present potential correlations between the measurements for a subset of the UV emission lines and the Lyα luminosity and the spectral slope β. These correlations are shown in Figures B.3, B.4, and B.5. Lastly, Figure B.6 shows the tentative relation between R_e and Δv_Lyα described at the end of Section 8.

Fig. B.1. Correlations between the individual components of the UV emission line doublets as described in Section 5. The MUSE objects (large symbols) are shown together with the literature comparison sample described in Appendix C (small dots) and color coded according to redshift.

Fig. B.2. Similar to Figure 6, but showing correlations between EW0 estimates as described in Section 6 instead of emission line fluxes. The MUSE objects (large symbols) are shown together with the literature comparison sample described in Appendix C (small dots) and color coded according to redshift.

Fig. B.3. Potential correlations between UV emission line flux of C III, He II, N V, and O III and log(LLyα) described in Section 7.

Fig. B.4. Potential correlations between EW0(C III), EW0(He II), EW0(N V), and EW0(O III) and log(LLyα) described in Section 7.

Fig. B.5. Potential correlations described in Section 7 between the O III emission line flux (left) and EW0(O III) (right) and the spectral slope β.

Fig. B.6. Tentative correlation described in Section 8 between the Lyα velocity offset and the effective radius, Re, of the LAEs.

Appendix C: Literature collection of UV emission line fluxes

In several of the plots in this paper the estimated UV emission line fluxes, EWs, line flux ratios etc. from the FELIS template matches to the MUSE-Wide, MUSE UDF mosaic and MUSE UDF10 spectra are compared to measurements collected from the literature. These UV emission line measurements were collected from Amorín et al. (2017, amo17), Bayliss et al. (2014, bay14), Berg et al. (2016, 2018, 2019a, 2019b, ber19), Christensen et al. (2012, chr12), Ding et al. (2017, din17), Du et al. (2020, du20), Erb et al. (2010, erb10), Herenz et al. (2020, her20), Hutchison et al. (2019); Jung et al. (2019, hut19), Jiang et al. (2020, jia20), Laporte et al. (2017, lap17), Le Fevre et al. (2019, lef19), Mainali et al. (2018, mai18), Mainali et al. (2020, mai20), Malkan et al. (1996, mal96), Marques-Chaves et al. (2020, mar19), Matthee et al. (2017, mat17), Nanayakkara et al. (2019, nan19), Ravindranath et al. (2020, rav20), Richard et al. (2021, ric21), Rigby et al. (2014, rig14), Rigby et al. (2015, rig15), Saxena et al. (2020, sax20), Schmidt et al. (2016, sch16), Schmidt et al. (2017), Mainali et al. (2017, sch17), Schaerer et al. (2018); Izotov et al. (2018, sch18), Senchyna et al. (2017, sen17), Senchyna et al. (2019, sen19), Shapley et al. (2003, sha03), Shibuya et al. (2018, shi18), Smit et al. (2017, smi17), Stark et al. (2014, sta14), Stark et al. (2017, 2015b, 2015a, sta15) Tang et al. (2021a, tan21), Vanzella et al. (2016, 2017, 2020, van20), and Wofford et al. (2021, wof21). The line fluxes and EWs collected for this literature sample are provided with this paper in the catalog described in Table C.1. The short key in the “reference” column of this table is indicated after each of the literature references. The redshift distribution of this sample of objects is shown as the dotted black histogram in Figure 2 where it is compared to the MUSE samples studied in this paper. Figure C.1 shows four projections of the multidimensional rest-frame EW space for the data sets in the literature catalog including plots of EW₀(C III), EW₀(C IV) for LAEs (bottom panels).

Fig. C.1. Projections of the rest-frame EW space spanned by EW0(He II), EW0(O III), EW0(C III), EW0(C IV) and EW0(Lyα) from the collection of UV emission line measurements from the literature plotted throughout this paper. The bottom panels reproduce the corresponding panels from Figure 10 and the top right panel reproduces the top right panel of Figure B.2. All panels exclude the MUSE measurements presented in this paper, but assign different symbols to the literature data included in Table C.1 according to the legend at the top. The points are color coded according to redshift and 3σ limits are indicated by the gray arrows. The diagonal curves correspond to the 10:1, 3:1, 1:1, 1:3, and 1:10 relations.

Appendix D: The full PIM-PDF of the NEOGAL and BPASS models

Figure D.1 shows the full prior PIM-PDFs of the physical parameters sampled by the NEOGAL and BPASS-based photoionization models described in Section 11. It is this initial photoionization model parameter space (listed in Table 5) that is constrained by the observations and measurements from the TDOSE spectra of the MUSE objects studied in this work resulting in object PIM-PDFs similar to the ones shown in Figure 20.

Fig. D.1. Full distribution of the PhotoIonization Model Probability Density Functions (PIM-PDFs) of the NEOGAL star forming (blue), NEOGAL AGN (purple), BPASS single (yellow), and BPASS binary (green) photoionization models (we note that the BPASS models fully overlap in initial parameters sampled). It is these models that are used for the parameter inference based on the observational constraints from the UV emission line detections presented in Section 11. The actual values of the sampled individual parameters are provided in Table 5.

Appendix E: Plots of TDOSE spectra and FELIS UV emission line detections

Figures E.1–E.6 provide further examples of sources with FELIS UV emission line detections in the TDOSE spectra similar to Figure 3 representing the breadth of sources in the parent sample analyzed in this paper. For a detailed caption we refer to Figure 3. Several of the panels in these figures support the discussion of individual sources as noted in the respective captions.

Fig. E.1. Similar to Figure 3, but showing the MUSE-Wide GOODS-South AGN 104014050 (z = 3.66) (top) and 115003085 (z = 3.71, bottom, see also Norman et al. 2002).

Fig. E.2. Similar to Figure 3, but showing the UDF mosaic AGN 601381485 (z = 3.19) (top) and object 604992563 (z = 3.80) (bottom). The latter object has a lead line redshift blue-wards of the Lyα peak as indicated by the vertical gray lines. Even when assigning the objects redshift according to a fit of the Lyα profile, as described in Section 8, the offset of the Lyα emission with respect to systemic as traced by the C III still has an intriguing value of −145 km s−1.

Fig. E.3. Similar to Figure 3, but showing the objects 102014087 (z = 3.25) (top) and 601281436 (z = 3.09) (bottom) for which we show the PIM-PDFs in Figure 20.

Fig. E.4. Similar to Figure 3, but showing COSMOS MUSE-Wide objects 202021046 (z = 3.95) (top) and 221004004 (z = 2.15) (bottom). The latter showing continuum, He II emission and broad C IV emission characteristic of AGN.

Fig. E.5. Similar to Figure 3, but showing the two lower-z UDF mosaic objects 600921283 (z = 2.17) (top) and 605172634 (z = 2.62) (bottom) where all covered lines were clearly detected. Object 605172634 is the outlier described in Section 9.

Fig. E.6. Similar to Figure 3, but showing the two UDF10 LAEs 720500425 (z = 3.28) (top) and 721250731 (z = 3.72) (bottom).

All Tables

All Figures

Fig. 1. Location of the analyzed MUSE pointings in the GOODS-South (top), HUDF (center) and COSMOS (bottom) regions. The MUSE-Wide pointings are shown in orange, the UDF mosaic in blue and the UDF10 pointing in red. The location of the 2052 sources searched for rest-frame UV emission lines are marked with green circles. The gray scale images are the 3D-HST combined near-infrared detection images from https://3dhst.research.yale.edu/ (top and bottom) and the HLF F775W version 2.0 image from https://archive.stsci.edu/prepds/hlf/ (center).

In the text

Fig. 2. Redshift distributions of the sources studied here (thick green histogram), the subsample of them with UV emission lines detected by FELIS (filled dark green histogram), the comparison sample of UV line emitters from the literature described in Appendix C (dotted black histogram), the MUSE-Wide DR1 catalog by Urrutia et al. (2019, thin orange histogram), and the independent MUSE-Deep UDF mosaic and UDF10 catalog presented by Inami et al. (2017, thin blue histogram). The MUSE redshift desert where neither [O II] λ3726 + [O II] λ3729 (0.3 < z < 1.5) nor Lyα (z > 2.9) emission is available for redshift identification in the MUSE wavelength range is marked by the gray band.

In the text

Fig. 3. Examples of TDOSE extracted 1D MUSE spectra (blue) of the two LAEs 201004081 (z = 3.08; top) and 602922055 (z = 3.44; bottom) with their 1σ error spectrum indicated by the blue-shaded region. Further examples can be found in Appendix E. Each of the smaller panels shows a UV emission line region as marked by the black boxes in the bottom panels, which show the full MUSE wavelength range. Vertical gray lines mark the approximate location of emission features based on the redshift of the LSDCat lead line. No continuum subtraction has been performed on the spectra. The UV emission lines detected by the software FELIS (Sect. 5) are marked by the red curves. The solid red curves show FELIS detections with S/N > 3 and a velocity offset below 1000 km s−1 deemed trustworthy in the visual inspection process. Dashed red curves show less secure detections, which were however deemed reliable. Dotted red curves show S/N > 3 and |Δv|< 1000 km s−1 FELIS detections classified as unreliable. In each panel 2× the median data cube noise in the respective fields (MUSE-Wide, MUSE mosaic and UDF10) are shown for reference as the black filled curves offset to the bottom of each panel for clarity. These median noise spectra are also shown in Fig. A.2.

In the text

Fig. 4. Percentage of objects with detected secondary UV emission lines red-wards of Lyα for the full sample (circles) and the subsample of LAEs (squares). The percentages are provided for all fields studied, the MUSE-Wide data (1 h depth), the UDF mosaic data (10 h depth) and the UDF10 pointing (31 h depth). The error bars indicate the 95% confidence Clopper–Pearson interval of the binomial distribution for k UV line detections given the n available objects searched (see Table 3).

In the text

Fig. 5. Distribution of rest-frame UV emission line fluxes (x-axis) detected at a FELIS S/N (y-axis) above 3σ (horizontal dashed line) in the MUSE-Wide (filled circles), the UDF mosaic (filled diamonds), and the UDF10 (filled x’s) data sets studied in this paper. These detections are compared to the collection of detections from the literature (small dots) described in Appendix C. Here S/N represents the emission line S/N quoted in the original references. All points are color coded according to their redshift. The solid histograms show the subset of the objects from this work. The dotted histograms show the distribution of measurements when including the literature data.

In the text

	Fig. 6. Similar to Fig. 5, but showing correlations between the > 3σ C III detections and the other secondary UV emission line fluxes from MUSE-Wide, the UDF mosaic, and UDF10 shown as large symbols. The literature measurements described in Appendix C are shown by the small dots. The dashed lines show the one-to-one relations to guide the eye.
In the text

	Fig. 7. AB magnitude distributions of the measured HST magnitudes used to infer the underlying continuum flux density when estimating the EW0 of the UV emission lines targeted in this study.
In the text

	Fig. 8. Rest-frame equivalent width (EW0) of the UV emission lines detected by FELIS as a function of the continuum magnitude used to estimate the flux density at the emission line location. As in Fig. 5, filled circles, diamonds and x’s correspond to objects from the MUSE-Wide, the UDF mosaic and the UDF10 samples, respectively. Each point is color coded according to each object’s redshift.
In the text

Fig. 9. Detected emission line fluxes as a function of rest-frame equivalent width of Lyα. The measurements presented here from MUSE-Wide, the UDF mosaic and the UDF10 are shown as large symbols. Measurements from the literature are shown by small dots (see Appendix C). All symbols are color coded according to source redshifts and are shown with 1σ error bars. The solid histograms show the subset of the objects from this work, whereas the dashed histograms include the estimates from the literature collection.

In the text

Fig. 10. UV emission line EW0 estimates of LAEs as a function of EW0(Lyα). The measurements from the current study are shown as large symbols (circles from MUSE-Wide, diamonds from the UDF mosaic, and x’s from UDF10). Measurements from the literature are shown by small dots (see Appendix C). All symbols are color coded according to source redshifts and are shown with 1σ error bars. The diagonal curves correspond to the relations 10:1, 3:1, 1:1, 1:3, and 1:10 for reference. Solid histograms show objects presented in this study, whereas dotted histograms include the literature samples.

In the text

Fig. 11. Velocity offsets for the individual FELIS template matches of the detected UV emission lines with respect to the object redshift based on the “lead line” of the MUSE LSDCat redshift determinations. For LAEs the Lyα redshift based on the line profile fits by Kerutt et al. (2021). From top left to bottom right the velocity offset for C IV, He II, O III, Si III, and C III with respect to the lead line are shown. For objects with z ≳ 2.9 (vertical dashed line) the lead line of essentially all objects is Lyα. The color coding shows the significance of the FELIS detection in terms of the signal to noise ratio, S/N(FELIS). All Δv estimates were assigned an uncertainty of ±100 km s−1 according to the precision uncertainty coming from the MUSE spectral resolution according to the discussion in Appendix A.3.

In the text

Fig. 12. Velocity offsets of Lyα determined with respect to the C III detections as a function of the Lyα rest-frame equivalent width estimated for a fixed UV slope of β = −1.97. The gray dot-dashed and dotted lines show the empirical relationship between ΔvLyα and equivalent width for z ≈ 3 objects presented by Nakajima et al. (2018) and Adelberger et al. (2003), respectively. The individual points are color coded according to redshift.

In the text

Fig. 13. Velocity offsets of Lyα as a function of the absolute UV magnitude for the MUSE LAEs with C III detections studied here shown by the large symbols and a collection of measurements at redshifts from z ≈ 2 to z = 7.73 collected from Kulas et al. (2012), Schenker et al. (2013), Erb et al. (2014), Stark et al. (2015a, 2017), Henry et al. (2015), Willott et al. (2015), Pentericci et al. (2016), Inoue et al. (2016), Bradač et al. (2017), Mainali et al. (2017), Verhamme et al. (2017), Sobral et al. (2018), Hutchison et al. (2019), Cassata et al. (2020), and Matthee et al. (2020a,b), shown by the small dots. To prevent cluttering the figure the estimated uncertainties are not shown. The median errors on the magnitudes and ΔvLyα are 0.2 mag and 50 km s−1. The dashed curves present the predicted median ΔvLyα as a function of UV magnitude for a range of redshifts (according to the color coding) from Mason et al. (2018b, corresponding to the black curve in their Fig. 2 from their Eq. (2)). This is based on an abundance matching model relating the absolute UV magnitude to the galaxy halo mass (Mason et al. 2015). Points and curves are color-coded according to object redshift and the color gradient highlights the importance of halo mass, that is available gas reservoir, in producing and determining ΔvLyα. The solid histograms show the distribution of objects from our study, whereas the dotted histograms include the sample of measurements from the literature.

In the text

	Fig. 14. Comparison of the velocity offset of Lyα with respect to the C III doublet detections and the Lyα velocity offset estimated based on the peak separation for double peaked LAEs and the Lyα FWHM for single peaked LAEs following the empirical relations presented by Verhamme et al. (2018). The dashed line indicates the one-to-one relation.
In the text

Fig. 15. Comparison of the velocity offset of Lyα and half the Lyα peak separation for objects with Lyα line profiles containing both a red peak and a blue bump (top) and the Lyα FWHM (bottom). These are the quantities used to predict ΔvLyα (Verhamme et al. 2018 approx.) on the x-axis in Fig. 14. The gray bands show the empirical relations presented by Verhamme et al. (2018). The small dots present a collection of literature measurements from Kulas et al. (2012), Leitet et al. (2013), Oesch et al. (2015), Hashimoto et al. (2015), Henry et al. (2015), Guaita et al. (2016), Izotov et al. (2016a,b), Stark et al. (2017), Verhamme et al. (2017, 2018), Yang et al. (2017), Matthee et al. (2020b), Muzahid et al. (2020), and Matthee et al. (2021). All points are color coded according to object redshift. The solid histograms show the measurements presented in this study, whereas the dotted histograms include the measurements from the literature.

In the text

Fig. 16. Electron densities (ne) estimated based on the C III (top) and Si III (bottom) doublet line flux ratios for a fixed electron temperature of Te = 104 K. The error bars (and limits) correspond to the 1σ uncertainties on the doublet flux ratios. For comparison, the theoretical curves (e.g., Osterbrock & Ferland 2006; Kewley et al. 2019a,b) for electron temperatures of 5000 K, 10 000 K and 20 000 K are shown as the solid red curves. The filled histograms count the number of objects where ne could be determined from the observed flux ratio. The open histograms add the subset of objects where the ne estimates are upper limits. The single open symbol to the right in the top panel indicates the C III flux ratio of object 605172634 (Fig. E.3) which results in an unconstrained electron density estimate for Te = 104 K.

In the text

	Fig. 17. Comparison of electron density estimates based on the C III (x-axis) and Si III (y-axis) from Fig. 16 for the 11 objects (all from the UDF mosaic and UDF10) with measurements from both doublets. The shown electron densities were estimated for a fixed electron temperature of Te = 104 K. The errorbars (and limits) corresponds to the 1σ uncertainties on the doublet flux ratios.
In the text

Fig. 18. Gas-phase abundances estimated using the third order polynomial fitting functions for the Si3-O3C3 (top panel) and He2-03C3 (bottom panel) diagnostics from Byler et al. (2020) as a function of source redshift. The large symbols are MUSE sources from the UDF and small dots are from the literature collection described in Appendix C. The solid histograms show the subset of the objects from this work, whereas the dotted histograms also include the sources from the literature. All points are color coded according to EW0(C III). The gray bands mark regions of super-solar gas-phase abundances at 12 + log10(O/H) > 8.69.

In the text

Fig. 19. Comparison of the gas-phase abundance estimates from Fig. 18 for objects with all relevant emission features detected. The large symbols are MUSE sources from the UDF and small dots are objects from the literature collection described in Appendix C. The solid histograms show the subset of objects from this work and the dotted histograms include the sample from the literature. The gray region marks values of super-solar abundances (12 + log10(O/H) > 8.69). Points are color-coded according to object redshift.

In the text

Fig. 20. Examples of PhotoIonization Model Probability Density Functions (PIM-PDFs) for the objects 102014087 (top) and 601281436 (bottom). The histograms in the different panels show the distribution of AGN (purple; Feltre et al. 2016), star forming (blue; Gutkin et al. 2016), BPASS binary (green) and single star (yellow; Xiao et al. 2018) photoionization model parameters satisfying the observational emission line flux ratio constraints listed in the left margin. The spacing of the parameter histograms’ bars indicate the discreteness of the sampling of each parameter space sampled by the models as listed in Table 5. The spectra of the two objects are shown in Fig. E.3. The limits on the observed flux ratios are all 3σ. The vertical dotted lines show the median values and the 68% confidence intervals of the distributions. The vertical dashed lines mark the mean values. For comparison the full unconstrained set of PIM-PDFs are shown in Fig. D.1.

In the text

Fig. 21. Mean value (x-axis) and normalized standard deviation (y-axis) of the ionization parameter (log10(U)) and the metallicity (Z) PIM-PDFs for the MUSE-Wide, the UDF mosaic, and the UDF10 objects studied in this work. The blue points represent the distribution of star forming models (Gutkin et al. 2016) matching the observations, the purple points indicate AGN models (Feltre et al. 2016), the green points show models based on BPASS binaries (Xiao et al. 2018), and the yellow points represent models of BPASS single stars (Xiao et al. 2018). The solid histograms show the distribution of the individual parameters. The vertical dashed line in the right panel marks the solar metallicity Z⊙ = 0.01524.

In the text

Fig. A.1. Overview of the FELIS template fit to a simulated spectrum of a C III emitter with a doublet flux ratio of F([C III])/F(C III]) = 1.4 corresponding to an electron density < 104 cm−3 (Osterbrock & Ferland 2006). The top panel shows the simulated spectrum (black) with the best-fit scaled template plotted on top (green). The second panel shows the S/N spectrum, that is flux/. The third panel from the top shows the template cross-correlation S/N calculated with Equation A.6. The maximum value is indicated by the vertical dashed red line in all four panels. The bottom panel shows the χ2 estimate of the template cross-correlation.

In the text

Fig. A.2. Median noise spectrum of two random MUSE-Wide pointings (green) and 70 GOODS-South and COSMOS MUSE-Wide pointings (blue – almost identical to the green curve). The latter spectrum was used to draw noisy flux values for the MUSE-Wide mock spectra that FELIS was tested on. For comparison, the median noise spectra of the deeper MUSE UDF mosaic and the UDF10 pointing are shown in red and gray. The median noise spectra shown in this panel are displayed by the black filled regions at the bottom of each panel in the spectral overviews in Figure 3 and Figures E.1–E.6.

In the text

Fig. A.3. Noise-free (top) and noisy (bottom) template matches to the mock spectrum shown in Figure A.1. The top panel of Figure A.1 shows the same as the bottom panel of this figure. The red curve shows the FELIS template match to the part of the spectrum in black (within the vertical dotted lines). The green curve shows the template match before applying the flux scaling. The part of the mock spectrum outside the considered wavelength range is marked in gray.

In the text

	Fig. B.1. Correlations between the individual components of the UV emission line doublets as described in Section 5. The MUSE objects (large symbols) are shown together with the literature comparison sample described in Appendix C (small dots) and color coded according to redshift.
In the text

	Fig. B.2. Similar to Figure 6, but showing correlations between EW0 estimates as described in Section 6 instead of emission line fluxes. The MUSE objects (large symbols) are shown together with the literature comparison sample described in Appendix C (small dots) and color coded according to redshift.
In the text

	Fig. B.3. Potential correlations between UV emission line flux of C III, He II, N V, and O III and log(LLyα) described in Section 7.
In the text

	Fig. B.4. Potential correlations between EW0(C III), EW0(He II), EW0(N V), and EW0(O III) and log(LLyα) described in Section 7.
In the text

	Fig. B.5. Potential correlations described in Section 7 between the O III emission line flux (left) and EW0(O III) (right) and the spectral slope β.
In the text

	Fig. B.6. Tentative correlation described in Section 8 between the Lyα velocity offset and the effective radius, Re, of the LAEs.
In the text

Fig. C.1. Projections of the rest-frame EW space spanned by EW0(He II), EW0(O III), EW0(C III), EW0(C IV) and EW0(Lyα) from the collection of UV emission line measurements from the literature plotted throughout this paper. The bottom panels reproduce the corresponding panels from Figure 10 and the top right panel reproduces the top right panel of Figure B.2. All panels exclude the MUSE measurements presented in this paper, but assign different symbols to the literature data included in Table C.1 according to the legend at the top. The points are color coded according to redshift and 3σ limits are indicated by the gray arrows. The diagonal curves correspond to the 10:1, 3:1, 1:1, 1:3, and 1:10 relations.

In the text

Fig. D.1. Full distribution of the PhotoIonization Model Probability Density Functions (PIM-PDFs) of the NEOGAL star forming (blue), NEOGAL AGN (purple), BPASS single (yellow), and BPASS binary (green) photoionization models (we note that the BPASS models fully overlap in initial parameters sampled). It is these models that are used for the parameter inference based on the observational constraints from the UV emission line detections presented in Section 11. The actual values of the sampled individual parameters are provided in Table 5.

In the text

	Fig. E.1. Similar to Figure 3, but showing the MUSE-Wide GOODS-South AGN 104014050 (z = 3.66) (top) and 115003085 (z = 3.71, bottom, see also Norman et al. 2002).
In the text

Fig. E.2. Similar to Figure 3, but showing the UDF mosaic AGN 601381485 (z = 3.19) (top) and object 604992563 (z = 3.80) (bottom). The latter object has a lead line redshift blue-wards of the Lyα peak as indicated by the vertical gray lines. Even when assigning the objects redshift according to a fit of the Lyα profile, as described in Section 8, the offset of the Lyα emission with respect to systemic as traced by the C III still has an intriguing value of −145 km s−1.

In the text

	Fig. E.3. Similar to Figure 3, but showing the objects 102014087 (z = 3.25) (top) and 601281436 (z = 3.09) (bottom) for which we show the PIM-PDFs in Figure 20.
In the text

	Fig. E.4. Similar to Figure 3, but showing COSMOS MUSE-Wide objects 202021046 (z = 3.95) (top) and 221004004 (z = 2.15) (bottom). The latter showing continuum, He II emission and broad C IV emission characteristic of AGN.
In the text

	Fig. E.5. Similar to Figure 3, but showing the two lower-z UDF mosaic objects 600921283 (z = 2.17) (top) and 605172634 (z = 2.62) (bottom) where all covered lines were clearly detected. Object 605172634 is the outlier described in Section 9.
In the text

	Fig. E.6. Similar to Figure 3, but showing the two UDF10 LAEs 720500425 (z = 3.28) (top) and 721250731 (z = 3.72) (bottom).
In the text